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ABSTRACT 



Three different approaches to the problem of imple- 
menting a reduced-order, sub-optimal Kalman filter for a 
discrete, linear stochastic process, with time-invariant 
dynamics, are presented. 

A first method. A, is based upon the partitioning of 
the system dynamics. A second method, B, is implemented 
using matrix pseudo-inversion and a third method, C, is 
based upon reduction of the original process to one of 
lower order using the dominant roots of the system. An 
expression for the performance degradation in method A is 
derived. In method B, expressions for the sub-optimal 
estimation error, and sub-optimal variance of estimation 
error are derived. 

The several methods are applied to a fourth-order 
process for illustration. 



JAVAL POSTuKrtDJAfE SCHOOL 
.lONTEREY, CALIF. 93940 



TABLE OF CONTENTS 

I. INTRODUCTION 7 

II. THE DISCRETE (OR DISCRETIZED) LINEAR iq 

I 

STOCHASTIC PROCESS WITH TIME- INVARIANT 
DYNAMICS 

III. IMPLEMENTATION OF A REDUCED DISCRETE, 16 

LINEAR KALMAN FILTER BY PARTITIONING 

THE DYNAMICS OF THE ORIGINAL PROCESS 

A. DERIVATION OF THE REDUCED FILTER 16 

EQUATIONS 

B. REDUCED FILTER DEGRADATION 24 

C. EXAMPLE 1 27 

IV. IMPLEMENTATION OF A SUB-OPTIMAL LINEAR 43 

FILTER FOR DISCRETE DYNAMIC PROCESSES 

USING MATRIX PSEUDO-INVERSION 

A. DERIVATION OF THE SUB-OPTIMAL FILTER 44 

EQUATIONS 

B. SUB-OPTIMAL ESTIMATES, ESTIMATION 54 

ERROR AND DERIVATION OF A LINEAR 

MATRIX DIFFERENCE EQUATION TO COMPUTE 
THE COVARIANCE OF SUB-OPTIMAL 
ESTIMATION ERROR 

C. EXAMPLE 2 59 

V. IMPLEMENTATION OF A REDUCED DISCRETE 81 

KALMAN FILTER BY USING THE DOMINANT ROOTS 
OF THE ORIGINAL PROCESS 



3 



A. THE DOMINANT ROOTS APPROACH 81 

B. EXAMPLE 3 92 

VI. DISCUSSION OF RESULTS 104 

APPENDIX A. MATRIX PSEUDO-INVERSION 106 

APPENDIX B. COMPUTER PROGRAM FOR PARTITIONING 109 

AND DOMINANT ROOTS METHODS 

APPENDIX C. COMPUTER PROGRAM FOR PSEUDO- 119 

INVERSION METHOD 

LIST OF REFERENCES 130 

INITIAL DISTRIBUTION LIST 131 

FORM DD 1473 133 



4 



ACKNOWLEDGEMENT 



The author wishes to express his deepest gratitude to 
Dr. James S. Deraetry for his most valuable guidance and 
patience in this work. He also wishes to thank his wife, 
Aida, for her constant encouragement and understanding. 



5 



I. INTRODUCTION 

In any modern control problem involving filtering of 
discrete or continuous stochastic linear ■ processes , it is 
clear that the computational requirements associated with 
the implementation of a Kalman filter increase consider- 
ably as the dimension of the system increases. A simpler 
(reduced) configuration of the filter with sub-optimal 
performance may be acceptable, mainly if some preselected 
r states (where rcn, and n represents the order of the 
original system) are of primary interest in the estimation 
process and if the performance degradation incurred when 
using the reduced model is not severe. 

With this idea of implementing a reduced, discrete, 
sub-optimal filter, based upon the optimal Kalman estima- 
tion theory (Ref. 1.),. a library search was done, yielding 
two papers of particular interest which present different 
approaches to the problem of reduced estimation for con- 
tinuous dynamic processes. Meditch (Ref. 3) in 1964, de- 
veloped a theory of sub-optimal filtering for continuous 
dynamic processes, in which the process vector is divided 
by means of matrix pseudo inversion, with separate filters 
designed for each part. The sub-optimal estimates are ob- 
tained by combining the several filter outputs. As 
Meditch points out in his paper, an important feature of 
his continuous formulation is that it provides a powerful 
tool for preliminary analysis of sub-optimal filters. 
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Huddle and Wismer (Ref. 4) in 1968 based a sub-optimal 
filter design on a "primary system", which is implemented 
by dividing the original process into two separate systems, 
using partitioning of the continuous dynamics of the sys- 
tem, and neglecting the remaining "secondary system". 

Another successful technique for simplifying filter 
computations has been investigated by Aoki and Huddle 
(Ref. 5), wherein a procedure is developed to estimate the 
state vector of a discrete stochastic system, when con- 
straints are imposed on the number of memory elements of 
the estimator. 

This thesis presents both an extension of two of the 
methods referred to above (Ref. 3 and Ref. 4) to the discrete 
case, and a third approach to the same problem based upon 
the use of the dominant roots of the full process. 

The study of the three different approaches to the 
reduced sub-optimal filter problem is presented in a co- 
herent form as follows. Chapter II describes the general 
discrete linear stochastic process. In Chapter III, the 
work of Huddle and Wismer (Ref. 4) is extended to the dis- 
crete case, and the performance degradation of the re- 
duced filter is obtained. An example using a fourth- 
order process is solved estimating only two states of the 
full system, in order to illustrate the method. In Chap- 
ter IV, the Meditch (Ref. 3) reduced filter is adapted to 
the discrete case, and recursive equations for the sub- 
optimal estimation error and sub-optimal covariance of 
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estimation error are derived and used as a means of com- 
paring the sub-optimal with optimal performances. The 
same fourth-order system used in Chapter III, is used here 
as an example, and again, only two states are estimated. 

In Chapter V, another reduced filter is implemented, using 
the dominant roots of the complete system in order to pro- 
duce the dynamics of the reduced filter. Once more, the 
technique is illustrated using the original plant of the 
previous examples. 

In the three examples presented as illustrations, the 
same complete process is used in order to compare the re- 
sults of the three different methods. In each one of the 
examples, a Monte Carlo simulation of the full process was 
performed, and optimal Kalman estimation of the whole 
state vector was accomplished in order to establish the 
performance degradation of the reduced filters. 

Appendix A presents the basic principles of matrix 
pseudo-inversion; more complete treatments of this sub- 
ject appear in Penrose (Ref. 8) and Werther (Ref. 9). 

Appendix B and C include the computer programs used 
in solving the examples. These programs were not intended 
to be computationally efficient either in time, or in 
storage. 
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1 1 . THE DISCRETE (DISCRETIZED) LINEAR 



STOCHASTIC PROCESS WITH 
TIME- INVARIANT DYNAMICS 

A linear, continuous stochastic dynamic process may 
be described in general (Ref. 7) by the vector matrix dif- 
ferential equation 

x(t) = A(t)x(t) + B(t)' u(t)+ w(t)j (2.1) 

with output 

y_(t) = H(t)x(t) (2.2) 

where: x(t) is a (n x 1) vector, representing the state 
variables of the system 

u(t) is a (m X 1) vector (m n), representing the 
deterministic input to the system 

w(t) is a (m X 1) vector representing the stochas- 
tic input to the system 

A(t) is a (n X n) system matrix 

B(t) is a (m X n) distribution matrix 

y^(t) is a (p X 1) vector of system outputs 

H(t) is a (p X n) measurement matrix. 

It will be assumed that measurement of the output 
can be achieved only with the addition of noise, so the 
measurement equation is defined as 
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z( t ) = ^( t ) + v( t ) , 



(2.3) 



where v(t) is a (p x 1) vector representing the stochastic 
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measurement noise and ^(t) is a (p x 1) vector of measured 



outputs. 

The process is visualized better by means of Fig. 2.1, 
where the double lines indicate a signal flow character- 
istic of all matrix block diagrams. 

It is assumed for simplicity that the process is sta- 
tionary (time-invariant); that is, all the elements of 
the A(t), B(t) and H(t) matrices are constants. It is 
further assumed that the process is to be sampled at dis- 
crete times and that u(t) and w(t) are constant vectors 
during each sampling period. The general solution for 
Eq. (2.1) under these conditions is expressed as 



If the sampling time is assumed to be T seconds, Eq. 




(2.4) may be rewritten as 




(2.5) 



o 



Defining 



^(T) = r-AT 



(2.6a) 
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Fig. 2.1 - Block diagram for a continuous stochastic 
process with time-varying dynamics. 




Fig. 2.2 - Block diagram for a discrete stochastic 
process with t ime- invariant dynamics. 
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and 



r (T) = B , (2.6b) 

o 

where $(T) (n x n) and P (T) (n x tn) are called the dis- 
crete state transition and distribution matrices, respect- 
ively. The discretized^ form of Eq. (2.1) becomes in the 
stationary case 



x(k+l) = J(T)x(k) + r(T)[u(k) + w(k)j . (2.7) 



^(T) and r (T) may be obtained by application of the 
Laplace transformation, i.e.. 



$(T) =o^"^[(sI 




(2.8a) 



1 



Distinction is made between a truly discrete system, and 
one which is continuously operating but driven and measured 
at discrete times. The linear difference equation for the 
truly discrete system will be 



x(k+l) = Aj^x(k) + u(k) + w(k) 



The terminology "discrete”, 
distinction from this point 



"discretized" is used without 
on . 
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(2.8b) 



r(T) 



_ p-l ^(s)B 



t = T 



The discretized, linear stochastic process with time- 
invariant coefficients will hereafter be represented by 
the model (see Fig. 2.2) 

x(k+l) = ^x(k) + r|^u(k) + w(k)j (2.9) 

with measurement equation 



^(k) = Hx(k) + v(k) . (2.10) 

Eqs. (2.9) and (2.10) will be used throughout this 
thesis as the model of the complete discrete system for 
which reduced sub-optimal filters are implemented. 
Furthermore, w(k) and v(k) will be considered to be vec- 
tors of independent, gaussian random processes (white 
noise) with zero means and fully described by the follow- 
ing statistics: 

E[w(k)| =0 for allk, (2.11) 



E,w(k)*w^(i) = E^w(k)j •E^w^(i)j =0 k i,(2.12) 



r T 1 ' ^ ‘ 

E|w(k)*w (k)' = known' quantity for all k , (2.13) 
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for all k , 



(2.14) 



Ejv(k)j = 0 

E[v(k)-v^(i)] = E[v(k)] •E[v'^(i)] =0 k 7 ^ i, (2.15) 

and 

E[v(k) • v^(k)| = R for all k , (2.16) 

where R is a (p x p) matrix of covariance of measurement 
noise. 
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III. IMPLEMENTATION OF A REDUCED DISCRETE, 



LINEAR KALMAN FILTER BY PARTITIONING 
THE DYNAMICS OF THE ORIGINAL PROCESS 

This chapter presents a discrete solution to the prob- 
lem of implementing a sub-optimal filter by means of par- 
titioning the original discrete system (Fig. 2.2) into two 
reduced models, a "primary system" model which will include 
only a selected subset of the state vector components whose 
estimation is of primary interest, and a "secondary sys- 
tem" which will be formed by the remaining state varibales. 
Then, the primary plant will have reduced dimensions and 
the implementation of a reduced Kalman filter based upon 
this primary model will be possible. 

A. DERIVATION OF THE REDUCED FILTER EQUATIONS 

It is assumed that the system dynamics are represent- 
ed by Eq . (2.9) and that the statistics of the forcing and 
measurement noise are fully described by Eqs. (2.11) 
through (2.16). For simplicity it will be further assumed 
that there is no deterministic forcing function u(k). The 
stochastic process may then be expressed as 

x(k+l) = $ 2 S(k) + Fw(k) , (3.1) 

with measurement 

^(k) = Hx(k) + v(k). 
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By matrix partitioning, the discrete dynamics of the 
full process may be rewritten as 



Xj^(k+1) 


rH 

t 


' ^ ll 




X 2 (k+1 ) 


:'lc2 

L 


1 $2 


1 

L 



Xi(k) 



q I 0 



(k )1 

^2 



X2(k) 



0 



Wo (k ) 



J 



( 3 . 2 ) 



With the full state vector being an (n x 1 ) matrix, 

then : 

represents an (n^^ x 1 ) state vector which includes 
all those states of particular interest in the estimation 
process ; 

X2 represents an (n2 x 1 ) vector of the remaining 
state vector components (n^^ + n2 = n); 

is an (n^^ x n^) matrix, representing the dynamics 

of Xi', 

$2 is an (n2 x n2) matrix, representing the dynamics 
of j<2 > 

and ^^2 ^ ^ ’^l^ matrices re- 

spectively, indicating the coupling between the primary 
and secondary systems; 



w^ is a 
functions ; 


(mj^ 


X 1 ) 


vector of 


forcing 


(white 


noise ) 


W2 is a 


(^2 


X 1 ) 


vector of 


forcing 


(white 


noise) 



functions with m^^ + m2 = m; and 

and 12 are (n^^ x (n2 x m2) distribution ma- 

trices respectively. 
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The observation matrix H is partitioned as follows, 



H = I H2j . 



The measurement equation may therefore be rewritten as 



It is assumed that all the measured outputs are in- 
cluded in the primary system, which makes H 2 = 0. This 
partitioning may not always be possible, and even if it 
is, it may require special rearrangement of the ^ vector. 

In Eq. ( 3 . 3 ) ^(k) is the original (p x 1) measured 
vector and is the (p x n^^) measurement matrix for the 
primary system. It may be seen from this special parti- 
tioning of the H matrix that the primary system vector 
x-^ will not only be formed by those states of primary in- 
terest in the estimation process, but will also include 
those states whose estimation may not be needed, but 
whose outputs are being measured. This will establish 
the condition that n^^^p. 

As a result of the preceeding discussion, Eq. (3.3) 
becomes 



z(k) = [h^ I H 2 ' f Xj^(k) 1 + v(k) . 

I I 



(3.3) 



I X2^^^) j 



z(k) = H^Xj^(k) + v(k) , 



(3.4) 
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which represents the measurement equation for the reduced 
system. 

The general system may be described as consisting of 
a primary system. 

xj^(k+l) = > (3.5a) 

and a secondary system described by, 

X2(k+1) = * (3.5b) 

Fig. (3.1) shows a complete representation of the parti- 
tioned system. 

In order to reduce computational requirements in the 
estimation of the preselected states, the design of the 
filter is based on the dynamics of the primary system, 
which can be represented as a simplified model of the re- 
duced process by 

x^(k+l) = r^Wj^(k) , ^ (3.6) 



1 

A bar over a quantity indicates a reduced -model quan- 
tity. 
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system. The deterministic component of the forcing 



with raeasurement 



z(k) = Hj^Xj^(k) + v(k) . 

In the simplified model of the primary system the term in- 
volving been neglected due to the fact that 

the secondary state vector considered at all 

in the implementation of the reduced filter. The conse- 
quences of using Eq. (3.6) instead of Eq . (2.5a), will be 
analyzed in Section B of this chapter. 

If the Kalman filter algorithm is applied using the 
dynamics of the full plant, the matrix equation for the 
gain G(k), the covariance of estimation error P(k/k) and 
the prediction covariance P(k+l/k) are given in the dis- 
crete case by 

G(k) = P(k/k-l)H’’^ H P(k/k-l)H'^ + r'’^ , (3.7a) 

P(k/k) = I - G(k)H P(k/k-l) , (3.7b) 



and 



P(k+l/k) = $ P(k/k)$^ + Q , 



(3.7c) 



where 



Q = r E j^w(k)w^(k) j r ^ . 

G(k) is an (n x p) matrix, P(k/k), P(k/k-l) and P(k+l/k) 
are (n x n)matrices. 
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The estimation equation is given by 



x(k/k) = ^x(k-l/k-l) + G(k)^^(k) - H^x(k-l/k-l) , 

(3.8) 



where x(k/k) is an (n x 1) vector of estimated states. 

By induction, the reduced sub-optimal expressions 
using the dynamics of Eq . (3.6) are 



Gi(k) = Pj^(k/k-l)HHHj^Pj^(k/k-l)H^ + r] 



-1 



(3.9a) 



P^(k/k) = - G^(k)H^'p^(k/k-l) , (3.9b) 



and similarly 



Pj^(k+l/k) = $^Pj^(k/k)$^ + , 



(3.9c) 



where 



Ql = qE[wj^(k)w^(k)]rJ' . 

p 

l(k/k) (n^ X n^) represents the sub-optimal covariance 
of estimation error and Pj^(k+l/k) (n^^ x n^^) is the sub- 
optimal prediction covariance of estimation error. 
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The sub-optirnal estimation equation is 



|j^(k/k) = $^X3^(k-l/k-l) + G3^(k)[z(k) - Hj^^j^Xj^Ck-l/k-l) . 

(3.10) 

Considering the problem of initialization for the 
discrete reduced covariance equations, it is seen as a 
good approach to use for P^(0/-1), the partitioned por- 
tion of P(0/-1) corresponding to Pj^j^(0/-1), or considering 
that 



P(0/-1) 



>11(0/-!) 


P^2(0/-l) 


P2i(0/-1) 1 

zx 1 


P22(0/-1) 



then 



Pj^CO/-!) = Pj^j^(0/-1) . 



For the filter initial states, 

Xj^(0/-1) = Xj^(0/-1) , 



where 



M V - 



x(0/-l) 



Xj^(0/-1) 

X2(0/-l) 
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Fig. 3.2 shows a complete block diagram of the reduced 
plant and filter. 



B. REDUCED EILTER DEGRADATION 

Because the dynamics of a simplified system repre- 
sented by 



were used for the implementation of the reduced filter, 
instead of the dynamics of the primary system 



the estimation of x-j^ is caused to be sub-optimal and in 
error. This error may be considered the performance de- 
gradation of the sub-optimal filter. 

Considering the reduced filter, its estimation error 
covariance is given by 



Xj^(k+1) = ^j^Xj^(k) + Fj^w^(k) 



(3.6) 



x^(k+l) = $]^X]^(k) + r^Wj^(k) , (3.5a) 



V(k) 




(3.11) 



where 



x^(k) = ^j^Xj^(k-l) + ^^j^X 2 (k-l) + rj^wj^(k-l) , (3.5a) 



and by proper use of Eq . (3.6) it may be rewritten as 



In the same way, 



x^(k/k) = $Jj^(k-l/k-l)+G^(k)[z(k) - H^$3^£^(k-l/k-l)] , 



(3.10) 



which by substitution of Eq . (3.4) for z(k), gives 
x^(k/k) = $^Xj^(k-l/k-l)+ Gj^(k) [H^Xj^(k) +v(k) -H^$j^Xj(k-l/k-l) . 



Further substitution of Eq. (3.12) for gives 



xj^(k/k) = ^3^£3^(k-l/k-l) + Gj^(k)Hj^ xj^(k) + 

- ^j^Xj^(k-l/k-l) + Gj^(k)v(k) . 



(3.13) 



Using Eqs. (3.12) and (3.13), the experimental estimation 
error covariance V(k) (Eq. 3.11) may be rewritten, after 
proper rearrangement, as 



r . r ^ 

V(k) = E| ,(I^ - Gj^(k)H^)(x^(k) - $j^X3^(k-l/k-l) 

2 '^ 

+ (l^-G3^(k)H^)$^^X2^^-l) + Gj^(k)v(k)] !. 

(3.14) 

By application of matrix algebra and elemental algebra of 
expectations, Eq . (3.14) may be expanded and after some 
manipulation gives 
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+ 
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V(k) = [l^-G^(k)Hjp^(k/k-l) + 



I^I^-G^(k)H^ 



$ciE[x2^k-l)] $ 






+E 



T 

i^(k/k-l)x^(k-l)] $Ji+f^iE[x 2 (k-l)eJVk-l)] [Ii-GiO<)hJ , 



(3.15) 



where 

e^(k/k-l) = Xj^(k) - Xj^(k-l/k-l) . 

The first term on the right hand side of Eq. (3.15) 
can be easily recognized as Pj^(k/k), the sub-optimal 
theoretical variance of estimation error based on the re- 
duced model. It is intuitively expected that P|^(k/k) > 
Pll(k/k). The remaining terms of Eq . (3.15) represent the 
further degradation of the reduced filter and are seen to 
be functions of that partitioned portion of $ which couples 
2<2 to the primary system, i.e., ^cl* s suggests that as 

the coupling between X 2 the primary system ($^.]^) de- 

creases, the quantity P^(k/k) more accurately represents 
the actual estimation error covariance. 

C . EXAMPLE 1 

The process used for this example was a fourth-order 
model, with continuous linear stationary dynamics subject- 
ed to a forcing function having both a deterministic and a 
random component, with the latter having zero mean and 
known variance. It was assumed that only one state vari- 
able was measured, subject to the presence of gaussian 
measurement noise with zero mean and known variance. It 
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was further assumed that both random sequences have sta- 
tistics as described by Eqs. (2.11) through (2.16). The 
reduced filter was implemented in order to estimate only 
two states (x-j^,X 2 ). 

The fourth-order continuous system is described by 





0 


1 


0 


0 




0 


x(t) = 


0 


0 


1 


0 


x(t) + 


0 




0 


0 


0 


1 




0 




-400 


-320 


-118 


-19_ 




400 



U( t )+w( t ) ; . 



A deterministic function u(t) =0.5 was used in the ex- 
ample. The random forcing function w(t) was assumed to 
have zero mean and variance equal to 0.01 square units, 
i.e., w(t) ; N(0.0, 0.01). 

The measurement equation is 

z(t) = 1 0 0 ojx(t) + v(t) , 



with the measurement noise v(t), considered to have zero 
mean and variance equal to 0.25 square units, i.e., 
v(t) ; N(0.0, 0.25). 

Before applying the reduced-filter technique dis- 
cussed in Section III, A, the continuous process was dis- 
cretized using the ideas presented in Chapter II, so the 
discretized system could be described as 

x(k+l) = $x(^^) r|u(k) + w(k) , 
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with measurement 



z(k) = [l 0 0 o] x(k) + v(k) , 

where 

(2.6a) 



was computed using the series expansion 



$ = I + AT + ^(AT) 




+ 



9 



and 



r 





B , 



(2.6b) 



was computed by means of 



I = (I + AT + 
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(AT) + 



)(IT - A ^ + 



a 2 if . 

^ 3.’ 



)B. 



A Monte Carlo simulation^ was performed using 200 
replications of the trajectory, each trajectory being 
formed by 50 samples with sampling time T = 0.05 seconds. 



1 

An IBM System/360-67 computer was used for solving all 
the examples presented in this thesis. A copy of the 
main computer program for this example is presented in 
Appendix B. 
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This simulation was carried out in order to produce approx- 
imate ensemble values for the experimental mean and experi- 
mental variance of estimation error in Xj^ and X 2 , as a func 
tion of time from t = 0, for both the complete, optimal fil 
ter and the reduced filter. 

For the Monte Carlo simulation each state was assumed 
to have uncorrelated gaussian initial conditions with zero 
mean and the following variances: 

2 

= 0.5 square units 

2 

if^X 2 = 2.0 square units 
2 

Cf^x^ = 2.0 square units 
2 

Cf^x^ = 4.0 square units. 

In order to obtain an unbiased estimate from the 
Kalman filter, the initial values assigned to the filter 
state variables were 



or 



x(0/-l) = e[x(0)] 



x(0/-l) 



0 

0 

0 

0 
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and the initial value of P(0/-l) was 



or 



P(0/-1) = e' |x(0) - x(0/-l)' [x(0) - x(0/-l)]^ j 



J 



r, 

) I 



= E^ l^x(O) - E x(0) j jx(0) - Ej^x(0)i 



= COV X(0) ; , 



0. 5 



P(0/-1) = 



2.0 



O 



o 



2.0 



4.0 



The reduced filter was implemented to estimate only 
two states (x^ and X 2 ), and for this purpose partitioning 
was applied as follows. The ^ and P matrices (from 
Table III-l) 



$(0.05) 



0.9990 0.0449 

-0.0066 0.9946 

-0.3644 -0.2981 

-12.2713 -10.1815 



' 0.0012 
1 


0.0000 


1 0.0480 

1 


0.0009 


1 0.8871 

1 


0.0307 


1 

i -3.9181 


0.3043 



r^(0.05) = [ 0.0001 0.0066 ! 0.3644 12.2713], 

I J 
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<1X1 



were partitioned as shown to obtain 






1 



0.9990 0.0449 

- 0.0066 0.9946 



and 



0.0001 

0 . 0066 



Correspondingly , 



Pl(0/-1) 



0.5 0 

0 2.0 



( 0 /- 1 ) = 



Hi = : 1 0_ 



and 



Qi 



0.000 0.000 

0.000 0.000 



(See Table III-l) 
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The reduced filter was then implemented using the set 
of Eqs. (3.9a,b,c) repeated here for convenience, 






V,: i.u 



G3^(k) = Pj^Ck/k^l^H^: H^P^(k/k-l)H^^ ^ ^ (3.9a) 



Pj^(k/k) = 'l - Gj^(k)H^jPj^(k/k-l) 



(3.9b) 



and 



P3^(k+l/k) 



= $3^P3^(k/k)$j^ + Q 



(3.9c) 



Attention should be called to the fact that the pre- 

X-* -I ' 

sence of a deterministic forcing function in this problem 
changes slightly the form of Ecj".’ (3.10), which can be re- 
written as 



Xj^(k/k) = $2^X|^(k-l/k-l)+ r^U2^(k7l)+ G^^(k) 



;(k) 



$l2ii(k-l/l<-l)+ r^uj^(k-l) 



to reflect the effect .of u(k) on the prediction step. 

The outputs of the computer simulation are shown in 
Tables III-2 through III-5, where nvmerical results versus 
time are tabulated for the gains, theoretical variance, 
experimental variance and experimental means of estimation 
error for both the complete and reduced filters. Figs. 
(3.3) through (3.6) show graphical representations com- 
paring complete and reduced filters performance. 
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Table III-l 



PHI MATRIX 

0.9999 0.0499 

-0.0066 0.9946 

-0.3644 -0.2981 

-12.2713 -10.1815 



0.0012 0.0000 

0.0480 0.0009 

0.8871 0.0307 

-3.9181 0.3043 



GAMMA MATRIX 

0.0001 

0.0066 

0.3644 

12.2713 



Q MATRIX 

0.0000 0.0000 0.0000 

0.0000 0.0000 0.0000 

0.0000 0.0000 0.0013 

0.0000 0.0008 0.0447 



0.0000 

0.0008 

0.0447 

1.5059 



VARIANCE OF FORCING NOISE - 0.01 



VARIANCE OF MEASUREMENT NOISE - 0.25 



INITIALIZATION 


- P(0/-1) 






0.5000 


0.0 


0.0 


0.0 


0.0 


2.0000 


0.0 


0.0 


0.0 


0.0 


2 . 0000 


0.0 


0.0 


0.0 


0.0 


4.0000 


INITIALIZATION 


- x(0/-l) 







0.0 



0.0 



0.0 



0.0 



T = 0.05 seconds. 
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Table III-2 



GAINS FOR FULL AND REDUCED FILTERS 



K 


G( If n 


G(2, 1 ) 


G( 3, 1) 


G(4,l» 


GR( If 1) 


GR(2.1) 


1 


0.6667 


0.0 


0.0 


0.0 


0. 6667 


0. 0 


2 


0.4071 


0.2333 


-0.2095 


-7.2845 


0.4071 


0. 2330 




0. 3097 


0.4075 


-0.6537 


-9.4633 


0.3102 


0.4259 


4 


0.2677 


0.5088 


-1.1719 


-9.0573 


0. 2708 


0. 5779 


5 


0.2477 


0.5314 


-1.6080 


-7.0288 


0. 2569 


0.6814 


6 


0.2349 


0.4601 


-1.8640 


-4.4015 


0.2540 


0.7347 


7 


0.2227 


0.4108 


-1.9266 


-1.9563 


0. 2542 


0. 7445 


8 


0.2087 


0.3175 


-1.8378 


-0.0721 


0. 2537 


0.7222 


o 


0. 1928 


0.2268 


-1.6557 


1.1851 


0.2512 


0.6798 


10 


0.1758 


0.1474 


-1.4305 


1.9077 


0. 2464 


0. 6274 


11 


0. 1584 


0.0825 


-1.1970 


2.2348 


0.2399 


0.5720 


12 


0.1412 


0.0322 


-0.9760 


2.2961 


0. 2321 


0.5177 


13 


0.1247 


-0.0051 


-0.7780 


2. 1931 


0. 2237 


0.4669 


14 


0. 1093 


-0.0315 


-0.6072 


1.9984 


0.2151 


0.4205 


15 


0.0950 


-0.0490 


-0.4637 


1. 7610 


0. 2064 


0. 3789 


16 


0.0819 


-0.0596 


-0.3457 


1.5120 


0. 1979 


0.341 8 


17 


0.0701 


-0.0647 


-0.2505 


1.2705 


0. 1897 


0.3090 


18 


0.0596 


-0.0660 


-0.1752 


1.0474 


0. 1819 


0.2800 


19 


0.0502 


-0. 0643 


-0.1167 


0.8480 


0. 1745 


0. 2544 


20 


0.0420 


-0. 0606 


-0.0723 


0.6746 


0. 1675 


0.2317 


21 


0. 0348 


-0.0556 


-0.0395 


0.5269 


0. 1609 


0.2115 


22 


0. 0286 


-0.0499 


-0.0161 


0.4036 


0. 1547 


0.1937 


23 


0.0233 


-0.0438 


-0.0001 


0.3026 


0. 1488 


0.1777 


24 


0.0188 


-0.0378 


0.0100 


0.2215 


0. 1433 


0. 1635 


25 


0.0151 


-0.0320 


0.0157 


0. 1576 


0. 1382 


0.1508 


26 


0.0120 


-0.0266 


0.0181 


0. 1084 


0. 1333 


0.1393 


27 


0. 0095 


-0.0217 


0.0181 


0.0714 


0.1287 


0.1290 


28 


0.0075 


-0.0174 


0.0166 


0. 0444 


0. 1243 


0.1197 


29 


0. 0059 


-0.0136 


0.0141 


0.0254 


0. 1203 


0. 1112 


30 


0. 0047 


-0.0105 


0.0111 


0.0127 


0. 1164 


0.1036 


31 


0.0037 


-0.0078 


0.0080 


0.0046 


0. 1127 


0. 0966 


32 


0.0031 


-0.0057 


0.0050 


-0.0000 


0.1092 


0.0902 


33 


0.0026 


-0.0041 


0.0023 


-0.0022 


0. 1059 


0.0844 


34 


0.0022 


-0.0028 


-0.0001 


-0.0028 


0. 1028 


0. 0790 


35 


0.0020 


-0. 0018 


-0.0022 


-0.0024 


0.0998 


0.0741 


36 


0.0018 


-0.0011 


-0.0038 


-0.0015 


0.0970 


0.0696 


37 


0. 0018 


-0.0006 


-0.0050 


-0.0004 


0. 0943 


0. 0654 


38 


0.0017 


-0.0003 


-0.0060 


0.0007 


0.0917 


0.0616 


39 


0.0017 


-0.0001 


-0.0066 


0.0017 


0. 0893 


0.0580 


40 


0.0017 


0.0000 


-0.0070 


0.0024 


0.0869 


0.0547 


41 


0.0017 


0.0001 


-0.0073 


0.0030 


0. 0847 


0.0517 


42 


0.0017 


0.0001 


-0.0074 


0.0033 


0. 0825 


0. 0488 


43 


0.0017 


0.0001 


-0.0074 


0.0034 


0.0804 


0.0461 


44 


0.0017 


0.0001 


-0.0073 


0.0033 


0.0785 


0. 0437 


45 


0.0017 


0.0000 


-0.0072 


0.0031 


0. 0765 


0.0413 


46 


0.0017 


0.0000 


-0.0071 


0.0028 


0.0747 


0.0392 


47 


0.0017 


-0.0000 


-0.0070 


0.0025 


0. 0730 


0.0371 


48 


0.0017 


-0. 0000 


-0.0068 


0.0022 


0.0713 


0.0352 


49 


0.0017 


-0.0001 


-0.0067 


0.0018 


0.0696 


0.0334 


50 


0.0017 


-0.0001 


-0.0066 


0.0015 


0.0680 


0. 0318 




G (4x1) 


represents 

filter. 


the matrix of gains 


for the 


optimal 




GR(2xl) 


represents 

filter. 


the matrix of gains 


for the 


reduced 
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Table III-3 



THEDRETICAL VARIANCE OF ESTIMATION ERROR 



K 


PK(1,1) 


PK(2t2 1 


PM3,3) 


PK(4,4) 


PR( 1, 1» 


PR( 2t2I 


1 


0.1667 


2.0000 


2.0000 


4.0000 


0. 1667 


2.0000 


2 


0.1018 


1.9603 


1.7605 


242.6279 


0. 1018 


1.9558 


3 


0.0774 


1. 8002 


2.3189 


346.3362 


0. 0776 


1.8684 


4 


0.0669 


1.5182 


3.9792 


261. 1479 


0. 0677 


1.7325 


5 


0. 0619 


1.1771 


5.7036 


142.2309 


0.0642 


1.5559 


6 


0.0587 


0. 8478 


6.6569 


60.6997 


0.0635 


1.3561 


7 


0.0557 


0.5745 


6.6825 


21.3209 


0. 0635 


1.1534 


8 


0.0522 


0.3709 


6.0419 


8.2223 


0.0634 


0.9640 


9 


0.0482 


0. 2307 


5.0782 


6.9855 


0. 062 8 


0.7970 


10 


0.0439 


0.1400 


4.0514 


9.3412 


0.0616 


0.6557 


11 


0.0396 


0.0847 


3. 1097 


11.6978 


0. 0600 


0.5391 


12 


0.0353 


0.0529 


2.3160 


12.9672 


0.0580 


0.4442 


13 


0. 0312 


0.0361 


1.6824 


13.1242 


0.0559 


0.3676 


14 


0.0273 


0.0282 


1. 1960 


12.4778 


0.0538 


0.3058 


15 


0.0237 


0.0255 


0.8337 


11.3650 


0.0516 


0.2559 


16 


0.0205 


0.0254 


0.5710 


10.0516 


0.0495 


0.2155 


17 


0.0175 


0.0264 


0.3855 


8.7171 


0. 0474 


0.182 7 


18 


0.0149 


0.0275 


0.2581 


7.4696 


0.0455 


0. 1557 


19 


0.0126 


0.0283 


0.1736 


6.3652 


0.0436 


0.1336 


20 


0.0105 


0.0286 


0.1199 


5.4260 


0. 0419 


0. 1152 


21 


0. 0087 


0.02 82 


0.0877 


4.6526 


0.0402 


0. 0999 


22 


0.0072 


0.0273 


0.0701 


4.0336 


0. 0387 


0.0870 


23 


0.0058 


0.0259 


0.0619 


3.5512 


0.0372 


0.0762 


24 


0. 0047 


0.0241 


0.0595 


3. 1854 


0.0358 


0.0670 


25 


0.0038 


0.0221 


0.0601 


2.9159 


0.0345 


0.0591 


26 


0.0030 


0.0199 


0.0619 


2.7237 


0.0333 


0.0524 


27 


0.002A 


0.0177 


0.0639 


2.5918 


0.0322 


0.0466 


28 


0.0019 


0.0155 


0.0654 


2.5057 


0.0311 


0.0416 


29 


0.0015 


0.0133 


0.0659 


2.4529 


0.0301 


0.0373 


30 


0.0012 


0.0114 


0.0654 


2.4236 


0.0291 


0.0335 


31 


0. 0009 


0.0096 


0.0640 


2.4099 


0.0282 


0.0302 


32 


0. 0008 


0.0081 


0.0618 


2.4059 


0.0273 


0.027? 


33 


0.0006 


0.0067 


0.0590 


2.4071 


0. 0265 


0.0247 


34 


0. 0006 


0.0056 


0.0557 


2.4105 


0.0257 


0.0224 


35 


0.0005 


0.0046 


0.0523 


2.4141 


0. 0250 


0.0204 


36 


0.0005 


0.0039 


0.0489 


2.4167 


0.0243 


0.0186 


37 


0. 0004 


0. 0033 


0.0455 


2.4177 


0.0236 


0.0170 


38 


0.0004 


0.0028 


0.0424 


2.4170 


0.0229 


0.0155 


39 


0. 0004 


0.0024 


0.0396 


2.4146 


0.0223 


0.0142 


40 


0.0004 


0.0022 


0.0371 


2.4109 


0.0217 


0.0131 


41 


0.0004 


0.0020 


0.0349 


2.4062 


0. 0212 


0.0120 


42 


0.0004 


0.0018 


0.0330 


2.4008 


0.0206 


0.0111 


43 


0.0004 


0.0017 


0. 0315 


2.3951 


0. 0201 


0.0102 


44 


0.0004 


0.0017 


0.0303 


2.3894 


0.0196 


0. 0094 


45 


0. 0004 


0.0017 


0.0293 


2.3839 


0.0191 


0.0087 


46 


0.0004 


0.0016 


0.0285 


2.3787 


0.0187 


0.0081 


47 


0.0004 


0.0016 


0.0279 


2.3741 


0.0182 


0.0075 


48 


0.0004 


0. 0016 


0.0275 


2.3700 


0.0178 


0.0070 


49 


0.0004 


0.0016 


0.0272 


2.3664 


0. 0174 


0.0065 


50 


0.0004 


0. 0016 


0.0270 


2.3635 


0.0170 


0.0060 



PK represents the diagonal elements of the theoretical 
covariance matrix for the optimal filter. 

PR represents the diagonal elements of the theoretical 
covariance matrix for the reduced filter. 
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Table III-4 



EXPERIMENTAL VARIANCE OF ESTIMATION ERROR 



K 


VK ( 1 , 1 ) 


VK ( 2 f 21 


VK ( 3 , 3 I 


VK ( 4 » 4 l 


VRI 1,1 I 


VP ( 2 , 2 I 


1 


0.1510 


1.9066 


1.9733 


4.0880 


0. 1486 


1.9041 


2 


0.1022 


1.8632 


1.7310 


225.4341 


0. 1008 


1.8617 


3 


0.0807 


1.7518 


2.2373 


326.4941 


0. 0794 


1.7496 


4 


0.0667 


1.5112 


3.8145 


257.3462 


0.0652 


1.5140 


5 


0.0677 


1 . 1743 


5.6299 


138.0874 


0.0663 


1.2011 


6 


0.0529 


0.8372 


6.0810 


60.2167 


0.0511 


0. 9078 


7 


0.0518 


0 . 5928 


6.5273 


23.2761 


0.0504 


0.7632 


8 


0.0487 


0.3939 


6.0208 


8.7271 


0.0492 


0.7093 


Q 


0.0464 


0.2513 


5.1600 


7.2510 


0.0497 


0.7732 


10 


0.0457 


0.1612 


4.4359 


9.6855 


0.0542 


0.9210 


11 


0.0420 


0.1002 


3.5490 


12.4134 


0. 0570 


1.0924 


12 


0.0391 


0.0617 


2.7683 


14.0120 


0.0559 


1. 1048 


13 


0.0368 


0.0401 


2.0579 


14.2776 


0.0615 


1.1834 


14 


0.0297 


0.0300 


1.3831 


14.0637 


0.0592 


1.2549 


15 


0.0264 


0.0256 


0.9971 


11.9529 


0.0587 


1.2175 


16 


0.0226 


0.0247 


0.6811 


12.3539 


0. 0593 


1.1740 


17 


0.0188 


0.0260 


0.4289 


9.7735 


0.0585 


1.0867 


18 


0.0153 


0.0261 


0.2742 


7.4905 


0.0616 


1.0225 


19 


0.0131 


0.0271 


0.1830 


6.0835 


0. 0639 


0.9195 


20 


0.0112 


0.0269 


0.1292 


5.5882 


0.0661 


0. 8078 


21 


0.0094 


0.0274 


0.0892 


5.2477 


0.0694 


0.7044 


22 


0.0074 


0.0252 


0.0707 


4 . 1038 


0. 0661 


0.5994 


23 


0.0061 


0.0240 


0.0651 


3.5547 


0.0661 


0.4966 


24 


0.0047 


0.0216 


0.0634 


2.8311 


0. 0649 


0,4195 


25 


0.0039 


0.0202 


0.0596 


3.2233 


0.0640 


0.3411 


26 


0.0032 


0 . 0185 


0.0569 


2.5204 


0.0593 


0.2732 


27 


0.0025 


0.0166 


0.0547 


2 . 7723 


0. 0555 


0.2257 


28 


0.0020 


0.0149 


0.0546 


2.8033 


0.0566 


0, 1949 


29 


0.0016 


0.0130 


0.0550 


2.6636 


0.0542 


0.1729 


30 


0.0013 


0.0114 


0.0581 


2.0651 


0. 0510 


0.1560 


31 


0 . 0010 


0.0098 


0.0584 


2.4213 


0.0474 


0. 1489 


32 


0.0008 


0.0083 


0.0567 


2.7746 


0.0459 


0.1498 


33 


0.0007 


0.0071 


0.0539 


2.5252 


0,0435 


0. 1523 


34 


0.0006 


0.0060 


0.0490 


2. 1855 


0.0391 


0.1563 


35 


0.0005 


0.0050 


0.0435 


2 . 5705 


0,0378 


0.1674 


36 


0.0005 


0.0042 


0.0420 


2.2536 


0.0363 


0. 1752 


37 


0.0004 


0.0035 


0.0437 


2. 5342 


0.0349 


0.1817 


38 


0.0004 


0.0029 


0.0446 


2.3070 


0.0350 


0. 1891 


39 


0.0004 


0.0025 


0.0421 


2.0512 


0,0378 


0.1947 


40 


0.0004 


0.0022 


0.0366 


2.4148 


0. 0393 


0.1985 


41 


0.0004 


0 . 0019 ' 


0.0300 


2.2901 


0.0403 


0,2 024 


42 


0.0004 


0.0018 


0.0256 


2.5146 


0.0405 


0.2035 


43 


0.0004 


0.0017 


0.0241 


2.2052 


0. 0411 


0.2022 


44 


0.0004 


0.0016 


0.0261 


2.3978 


0.0429 


0.2016 


45 


0.0004 


0.0015 


0.0300 


1.9640 


0.0439 


0.1990 


46 


0.0004 


0.0014 


0.0317 


1.9171 


0. 0476 


0.1957 


47 


0.0004 


0.0014 


0.0307 


2.2013 


0,0494 


0.1914 


48 


0.0004 


0.0014 


0.0299 


2.3214 


0,0500 


0.1835 


49 


0.0004 


0.0015 


0.0281 


2.2047 


0.0532 


0.1772 


50 


0 . 0004 


0.0016 


0.0229 


2.8952 


0.0570 


0.1698 



VK represents the diagonal elements of the experimental 
covariance matrix for the optimal process. 

VR represents the diagonal elements of the experimental 
covariance matrix for the reduced process. 
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Table III-5 



EXPERIMENTAL MEANS OF ESTIMATION ERROR 



K 


MEANEI 


MEANE2 


MEANE3 


MEANE4 


MEANEIR 


MEANE2R 


1 


0. 0496 


0.0500 


-0.0350 


-0.1543 


0.0496 


0.0500 


2 


0.0372 


0.0392 


-0.0572 


-0.6276 


0. 0372 


0.0392 


3 


0.0356 


0.0307 


-0.0903 


-0.7992 


0.0358 


0.0429 


4 


0.0380 


0.0271 


-0. 1340 


-0. 7676 


0.0391 


0.0662 


5 


0.0395 


0.0203 


-0.1637 


-0. 3713 


0.0425 


0.0973 


6 


0. 0333 


-0.0030 


-0.1219 


-0.0534 


0.0387 


0.1099 


7 


0.0299 


-0.0146 


-0.0909 


0.2345 


0.0387 


0.1346 


8 


0.0044 


-0.0565 


0.1352 


0.1478 


0.0122 


0.0810 


9 


-0.0017 


-0.0536 


0.1711 


0.0767 


0.0089 


0.1003 


10 


0.0049 


-0.0373 


0.1013 


0.1589 


0.0229 


0.1579 


ll 


0.0142 


-0.0265 


0.0168 


-0.0124 


0.0414 


0.2121 


12 


0.0021 


-0.0281 


0.0951 


-0.0444 


0.0257 


0.1751 


13 


-0.0028 


-0.0232 


0. 1136 


-0.1564 


0. 0209 


0.1623 


14 


-0.0072 


-0.0168 


0.1219 


-0.3115 


0.0155 


0. 1447 


15 


-0.0086 


-0.0108 


0.1078 


-0.4010 


0.0148 


0.1327 


16 


-0.0058 


-0.0083 


0.0729 


-0.3752 


0. 0232 


0. 1322 


17 


-0. 0028 


-0. 0080 


0.0511 


0.0015 


0.0320 


0.1267 


18 


0.0020 


-0.0112 


0.0313 


-0.0783 


0.0461 


0. 1255 


19 


0.0055 


-0.0150 


0.0153 


-0.0996 


0.0572 


0.1156 


20 


0.0044 


-0.0137 


0.0109 


-0.1052 


0.0511 


0.0790 


21 


0.0040 


-0.0136 


0.0142 


0.2029 


0. 0476 


0.0468 


22 


0.0014 


-0.0094 


0.0236 


0.1062 


0.0321 


0.0012 


23 


0.0022 


-0.0103 


0.0273 


0.0610 


0.0343 


-0.0198 


24 


0.0017 


-0.0089 


0.0262 


-0.0824 


0. 0284 


-0. 0492 


25 


0. 0011 


-0.0072 


0.0235 


-0.0285 


0.0199 


-0.0795 


26 


0.0008 


-0.0063 


0.0215 


-0.0504 


0. 0143 


-0. 1043 


27 


0.0007 


-0.0056 


0.0207 


-0.0033 


0.0100 


-0. 1258 


28 


0.0006 


-0.0048 


0.0263 


0. 1783 


0. 0050 


-0.1455 


29 


0.0001 


-0.0027 


0.0329 


0.1111 


-0. 0078 


-0.1700 


30 


0.0002 


-0.0014 


0.0292 


-0.2192 


-0.0114 


-0.1832 


31 


0.0003 


-0.0006 


0.0215 


-0.1151 


-0.0128 


-0.1929 


32 


0. 0004 


0.0002 


0.0170 


-0.0758 


-0.0172 


-0.2037 


33 


0.0005 


0.0008 


0.0080 


-0.2516 


-0.0234 


-0.2143 


34 


0.0004 


0.0011 


0.0012 


-0.0467 


-0.0344 


-0.2271 


35 


0.0004 


0.0012 


0.0022 


0. 0639 


-0.0455 


-0.2379 


36 


0.0005 


0.0014 


0.0016 


-0.0675 


-0.0514 


-0.2430 


37 


0.0006 


0.0014 


-0.0011 


-0. 0404 


-0.0574 


-0.2467 


38 


0.0006 


0.0012 


-0.0059 


-0.1685 


-0.0654 


-0.2506 


39 


0.0006 


0.0007 


-0.0093 


0.0404 


-0. 0711 


-0.2520 


40 


0.0006 


0.0003 


-0.0057 


0.0877 


-0.0785 


-0.2533 


41 


0. 0007 


0.0002 


-0.0006 


0.1207 


-0.0795 


-0.2494 


42 


0.0007 


0.0003 


0.0027 


0.0151 


-0. 0876 


-0.2489 


43 


0.0008 


0.0004 


0.0038 


0.0374 


-0.0882 


-0.2432 


44 


0.0008 


0.0006 


0.0026 


-0.0628 


-0. 0905 


-0.2381 


45 


0.0008 


0.0007 


0.0004 


-0. 0407 


-0.0978 


-0.2354 


46 


0. 0008 


0.0007 


0.0035 


0.1341 


-0.1020 


-0.2307 


47 


0.0007 


0.0010 


0.0075 


0.0162 


-0.1105 


-0.2277 


48 


0. 0008 


0.0014 


0.0069 


-0.0321 


-0. 1128 


-0.2212 


49 


0. 0009 


0. 0017 


0.0066 


0.0174 


-0.1143 


-0.2143 


50 


0.0010 


0.0020 


0.0053 


-0.0515 


-0. 1149 


-0.2071 



MEANEI through MEANE4 represent the mean of estiaation 
error in the complete state vector by means of optimal 
filtering. 

MEANEIR and MEANE2R represent the mean of estimation 
error in states and X 2 by means of a reduced filter. 
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Variance of estimation error 




Fig. 3.3 - Graphical plot of the variance of estima 
tion error in x-^ vs. time. Curves a and b show the 
variance of estimation error, both theoretical and 
experimental, for the optimal estimation process. 
Curves c and d show the corresponding theoretical 
and experimental variances for the reduced case. 
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Variance of estimation error 




Fig. 3.4 - Graphical plot of the variance of est. 
error in X2 vs. time. Curves a and b show the vari- 
ance of estimation error, both theoretical and ex- 
perimental, for the optimal estimation process. 
Curves c and d show the corresponding theoretical 
and experimental variances for the reduced case. 
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Mean of estimation error 



0.20 



0.15 

0.10 

0.05 

0.0 

-0.05 

- 0.10 

- 0.12 



'1 




1.0 1.5 2.0 

Time in seconds 



Fig. 3.5 - Graphical plot of the means of estimation 
error in vs. time. Curve a shows the mean of 
estimation error for the optimal Kalman filter, and 
curve b shows that for the reduced filter. 



41 



Mean of estimation error 



0.4 



0.3 




' ! 

I i 

0.0 0.5 1.0 1.5 2.0 2.5 



Time in seconds 



Fig. 3.6 - Graphical plot of the means of estimation 
error in X 2 vs. time. Curve a shows the mean of 
estimation error for the optimal Kalman filter, and 
curve b shows that for the reduced filter. 
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IV. IMPLEMENTATION OF SUB-OPTIMAL LINEAR 



FILTERS FOR DISCRETE DYNAMIC PROCESSES 
USING MATRIX PSEUDO-INVERSION 

In the previous chapter a complete discussion of the 
implementation of a discrete, reduced model Kalman filter 
was presented, where the idea of partitioning the dynamics 
of the original plant was principal. In this chapter, a 
different approach to the same problem is presented, 
based upon the work of Meditch (Ref. 3) for the continuous 
case. 

Applying the same idea of reduction of the original 
process, it is partitioned for this case into two or more 
reduced vector components (i = 1,2 ... M) , with each 
new vector containing the same number of states. This 
new partitioning is accomplished by using matrix pseudo- 
inversion (see Appendix A). Once the original system has 
been broken down into M^ reduced-order sub-optimal fil- 
ters, then the sub-optimal estimates for the original pro- 
cess are obtained by recombination of the outputs from all 
the sub-optimal filters. Section A presents the deriva- 
tion of the discrete sub-optimal filter formulation. 

These results are extended to Section B, where expressions 

1 

The order (n) of the original system is equal to M 
times the order (r) of the reduced plants, that is 
n = M • r. 
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for the sub-optimal estimates, theoretical sub-optimal 
estimation error and theoretical sub-optimal covariance of 
estimation error are derived. Finally in Section C an ex- 
ample is solved to illustrate the technique. 

A. DERIVATION OF THE SUB-OPTIMAL FILTER EQUATIONS 

As in Chapter III-A it is assumed that the system 
dynamics are represented by Eq. (2.9), without consider- 
ing the presence of the deterministic part of the forcing 
vector, that is 

x(k+l) = ^x(k) + Pw(k) , 

with measurement 

^(k) = Hx(k) + v(k) , 

where as before, w(k) and v(k) are random gaussian func- 
tions fully described by Eqs. (2.11) through (2.16). 

If, as discussed in Chapter I, it is desired to 
estimate only r states, where r<^n , then the original vec- 
tor x(k) can be partitioned or transformed into two or 
more vectors by using a set of constant matrices 

A^ , i = 1 , 2 . . . M 

with A^ having dimensions (r x n). Notice that each A^ 
must have r rows and n columns. Then it is possible to 
define a new set of reduced vectors 
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^^(k) = A^x(k) , 



i = 1,2 ... M 



(4.1) 



with Xi(^) representing a (r x l) vector. 

Now the concept of the Generalized Inverse^ of a 
matrix (pseudo-inversion) is applied to the implementa- 

4 - 2 

tion of the sub-optimal filters. A matrix A is con- 
sidered to be the pseudo-inverse of a rectangular matrix 
A if 

AA^A = A , (4.2) 

where in general Aisa(sxt) matrix and A"*” is a (t x s) 
matrix. 

The selection of the set of constant matrices A^ will 
have as a necessary requirement for the solution of the 
filter implementation and the proper recovery of the sub- 
optimal estimates from the combination of the filter out- 
puts, that 

M 

^A^A. = I , (4.3) 

i=l 



1 

See Deutsch (Ref. 10) and Penrose (Ref. 8) for a complete 
discussion about matrix pseudo-inversion. 

2 

A stands for the pseudo-inverse of A. 
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where I is the identity matrix. This requirement will 
be clarified after the derivation which follows. 

If Eq. (4.1) is premultiplied by A^, i = 1,2 ... M 
and summation over all the i's is performed, then 



M M 

^ 2^ AtA^x(k) 
i=l i=l 



from which it is clearly seen, applying the requirement 
established by Eq . (4.3), that 

M 

x(k) = ^ Aty^. (k) . (4.4) 

i=l 



This last equation provides the key for recovery of 
the original system state vector. 

Now if the original discrete system equation is pre- 
multiplied by A^, i = 1,2 ... M, 

A^2S(1<‘'"1) = A^^x(k) + A^Fw(k) , (4.5) 



and if Eqs. (4.1) and (4.4) are substituted in Eq. (4.5), 
the latter becomes 



M 

y^(k+l) = A^J ^ Aty^(k)+A^r w(k) , i = 1,2 ... M 

j=l 

(4.6) 
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Expression (4.6) represents a set of M reduced, dis- 
crete system equations, each of order r. 

The measurement process can be written assuming that 
all the measurements that were available for the complete 
filter are available too for the reduced filters. 

^(k) = Hx(k) + v(k) , 

becomes after replacement of x(k) by Eq . (4.4), 

M 

z(k) = H Y' At^.(k) + v(k) . (.47) 

j=l 

Fig. (4.1) will help to better visualize the partitioned 
system. 

In order to facilitate the following derivations, the 
optimal Kalman filter estimation Eq. (3.8) may be rewritten 
in a slightly different form, that is 

x(k/k) = $x(k-l/k-l) + G(k)|^z(k) - z(k/k-l)J , 

(4.8) 

where £(k/k-l) represents the expected value of ^(k) given 
all measurements up to time (k-1), or 

^(k/k-1) = E[^(k)j given ^(k-l) . 



47 




48 



tioned into M sub-systems by means of matrix pseudo- 



From Eq. (4.7), however, 



M 



z(k) H ^ At ■*■ X(l<) • 

j = l 



(4.7) 



Therefore , 



M 



z(k/k-l) = 



H At E! (k)l given z(k-l) 

‘ — ' J ^ J J 

j = l 



or 



M 



z(k/k-l) = 



H a| 2.(k/k-i) 

j=l 



(4.9) 



Considering further that 



x(k/k-l) = $x(k-l/k-l) , 



(4.10) 



where x(k/k-l) represents the expected value of x(k) given 
all measurements up to time (k-l), and premultiplying Eq. 
(4.10) by A^, 1,2 ... M, 



A^x(k/k-l) = A^$x(k-l/k-l) . 



(4.11) 



By extension of Eq. (4.1) 



Aix(k/k-l) = » 
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and similarly by proper use of Eq. (4.4) 



M 

x(k-iA-i) = ^ , 

j=l 



so Eq. (4.11) becomes 



M 

2^(k/k-l) = Aj 2,j(k-l/k-l) . (4.12) 

j = l 



Basing the implementation of the reduced filters upon 
the dynamics of the set of systems described by Eq. (4.6), 
a set of M reduced-filter estimation equations may be 
stated by induction from Eq. (4.8), as 



2^(k/k) = + G^(k) z(k) - z(k/k-l)| 



Replacing for 2j^(k/k-l) and 2(k/k-l) Eqs. (4.12) and (4.9), 
gives 



M 



M 



Y At£j (k-l/k-l)+G^(k) z(k)-H Y At2j(k/k-l) 



j=l 



j = l 



(4.13) 



with i = 1,2 ... M. 

t It. 

This expression represents a set of M r -order coupled 
reduced filters, whose outputs £^(k/k) will be recombined 
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in a proper way in order to get the sub-optimal estimates 
x(k/k), where x(k/k) is an (n x 1) vector (see Section 
IV-B). 



Referring to the optimal gain and theoretical vari- 
ance expressions for the Kalman filter (Eqs. 3.7a,b,c), 
the equivalent M subsets of sub-optimal equations for each 
one of the reduced filters may be expressed, by definition, 
as 



G. (k) = P. (k/k-l)HT [h.P. (k/k-l)H'!^+R 

-L i 1 L 1 1 



, -1 



(4. 14a) 



where 



= H (see Eq . 4.7) 



corresponds to the i^^ subsystem contribution to the 



measurement vector; 




(4. 14b) 



and 



P.(k+l/k) = $^P.(k/k)$T + Q. 



(4.14c) 



where 



$i = (see Eq. 4.6) 

= A^r E[w(k)w'^(k)] t'^aT 



= a.qaT 

1 
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and 



r. = A . r . 

1 1 

It is important to notice that although the M estima- 
tion equations and the measurement equation are coupled 

Cor all the M reduced filters, the gain and covariance 

t h 

equations for the i*" filter are truly uncoupled from the 
other j = 1,2 ... M filters, where j i. 

Eqs. (4.14) complete the implementation of the M sub- 
optimal filters. The filtering process will be visualized 
better by means of Fig. (4.2). 

The initialization values for the sub-optimal filters 
are derived from the statistics used to initialize the com- 
plete optimal filter in the following way. Given the 
values for P(0/-1) and knowing that 

P^(0/-1) = e[x.(0)xT(o)] 



and 



y_^(0) = A^x(O) , 

then 



p.(o/-i) = E :a.x(o)x'^(o)aT 



1 — ' 



IJ 



A.E 

1 



X 



(0)x’’(0)l aT 



A.P(0/-l)At . 
1 1 



(4.15) 



52 




53 



Fig. 4.2 - Block diagram representation of the M sub-optimal filters, showing 
the outputs from each reduced filter, (k/k). 



By the same reasoning, 



£.(0/-l) = A.x(0/-1) 



(4.16) 



will be used as the initialization values for the set of M 
reduced estimation Eqs. (4.13). 

B. SUB-OPTIMAL ESTIMATES, ESTIMATION ERROR AND SUB- 

OPTIMAL VARIANCE OF ESTIMATION ERROR 

In the preceding section, M reduced filters were im- 
plemented after breaking down the original process into M 
reduced (r x l) vectors. In order to obtain the sub-optim- 
al estimates of the original process, the outputs £^(k/k) 
from the reduced filters have to be recombined. This is 
easily accomplished by proper application of Eq. (4.4) to 
the filter's estimates, whereupon 



will represent the sub-optimal estimates of the original 



The remainder of this section is devoted to the de- 
velopment of a matrix difference equation for the itera- 
tive computation of the sub-optimal covariance of estima- 
tion error. 

If the error of sub-optimal estimation is defined as 



M 




(4.17) 



process, with x(k/k) representing an (n x 1) vector. 



e(k) = x(k) - x(k/k) 



(4.18a) 
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(4.18b) 



then using Eq . (4.17) it becomes 

M 

e(k) = 2s(k) - At ^.(k/k) , 

j=l 

or 

M 

¥(k+l) = x(k+l) ■ ^ • (4.18c) 

j=l 



Substitution of Eqs. (4.13) and (4.7) into Eq . (4.18c) 
gives 

M (- M 

e(k+l) = x(k+l) - X!! ^i^ 

j=l j=l 

M M 

+G^(k+1)[H ^ At (k+l)+v(k+l)-H^ ^ At2j(k/k)] 

j=l j=l 

M 

Replacing y Aty. (k+1) with x(k+l), and substituting 

J J ” 

j=l 

x(k+l) = $x(k) + Pw(k) from Eq. (3.1), then 



M 



e(k+l) = $x(^) r*w(k) - 



j=l 



M 

A.$ r A^2j(k/k) 

j=l 

M 



+G^(k+1) [H$ 2 s(k) + Hrw(k) + v(k+l) - ^ Aty;.(k/k)] 

j=l ^ ^ 
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Expansion of this expression for the sub-optimal error 
leads to 



M M 

e(k+l) = $x^^) + rw(k) - ^ At£j(k/k) 

j=l j=l 



M M 

- AtGj(k+l)H$x(k) - ^ AjGjCk+DHr w(k) 

j=l j=l 



M KM 

- ^ AtGj(k+l)v(k+l)+ ^AjGj(k+l)H$ ^At^jO^k). 

j=l j=l j=l 



Collecting terras, 



M 

e(k+l) = [$ - X] A|Gj(k+l)H^][x(k) 

j=l 



M 

^ A+^.(k/k)] 



j = l 



M 

+ |I- A^G. (k+l)H 



M 

rw(k)- y AtG^+lMk+l). 

“ ^ J J “ 

j = l 



Substitution of Eq . (4.18b) for e(k) gives 

M M 

e(k+l) = [l- , AtG^(k+l)H|$e(k)+[l - ^ AtGjCk+DnJr^^k) 

j=l ^ j=l 

M 

■ X AtG^(k+l)v(k+l). (4.19) 

j = i 
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Eq . (4.19) represents a theoretical form for the error 
of sub-optimal estimation which can be useful in establish- 
ing the steady-state estimation error for the sub-optimal 
filter. 

If the sub-optimal theoretical variance of estimation 
error is defined as 



^(k) = Ei¥(k)e'^(k)J 



(4.20) 



where ^(k) is an (n x n)matrix, then substitution of Eq . 
(4.19) for ¥(k) gives 



f 

( r 



S(k+1) = E 




M M 

^ A|Gj(k+l)H]$i(k)+ [l- ^ AtG.(k+l)H]rUk) 

j=l j=l ^ ^ 



M ^ T 

+ 1 ! 

- 2^ AjGj (k+l)v(k+l) same terms 

j=l ^ ^ 



Carrying out the indicated matrix operations term by term 
yields 

/ M M 

^(k+1) = E([l- ^ AtGj(k+l)H]$e(k)e^(k)$^[l-y’AjGj(k+l)H] 

I j=l j=l 



,1T 



M 



M 






I-) AtG.(k+l)Hlrw(k)w'^(k)r'^fl- y AtG.(kd-l)H, 

J J J— — J 

j=i j=i 

M M ^ 

+ )^A^Gj(k+l)v(k+l)v'^(k+l) ^GT(k+l)At'^+ pJoduds 
j=l j=l j 
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Given the statistics of the problem, all the cross products 
will yield zero expectation. S(k+1) then becomes 



^(k+1) = 



M ^ 

“ X A H ^ E [e(k)e'^(k)] [ I - V Aj Gj(k+1)H ] 

j = l ^ j^^l 



M M ^ 

+ [l - A^G-Ck+DH] rE[^k) w'^Ck)] r^[l- y A‘l'G.(k+l)H] 

M M 

+ '■ AtG.(k+l)E[v(k+l)v^(k+l)] ^ gT ( k+l)A't'^ , 

j=l ^ ^ j=l ^ ^ 



and finally 



M 

S(k+1) = \ l - Y AtG.(k+l)H $S(k)$^ 
^ J J 

j=l 



y aIg. ( k+i)H 1 

^ j j 

j=i 



M M 

+ \l - ' " AtG.(k+l)HjQ [l - "^AlG.(k+l)H 
*- J J - J J 

j=l j=l 



M 

+ ^ AtG^(k+l)R ' GT(k+l)At^ . (A. 21) 

j=i j=i 



Eq. (4.21) represents an iterative expression for the 
computation of the theoretical sub-optimal covariance of 
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estimation error. To initialize the computation of ^(k+1), 
the initial value of ^(0) should be chosen as 



^(0) = E e(0)e^(0) 

M 

= E x(0) - At2^(0/-1) 

j = l 

2 ^ 

= E ^2S(0) - 2S(0/-1) 



= P(0/-1) . 



2 



The discrete matrix expression for the theoretical 
sub-optimal covariance of estimation error "S^(k), may be 
used as a reference to evaluate the performance of the re- 
duced set of filters, i.e., how close to the "true" values 

are the sub-optimal estimates. The solution for S(k) 
should be compared with the experimental variance of sub- 

optimal estimates. 

C. EXAMPLE 2 

The model used to investigate and to illustrate the 
technique presented in this chapter, was exactly the same 
fourth-order plant used in Example 1, subjected to the 
same deterministic and random forcing functions and with 
precisely the same measurements. The continuous system 
was discretized using the same sampling time T = 0.05 sec.. 
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and therefore the same $ and P matrices were obtained, 
arriving at the discretized system described by 

x(k+l) = ^ > (2.7) 

with measurement 

^(k) = Hx(k) + v(k) . 

The goal of this example, as in Example 1, was to 
estimate only two states of the fourth-order system 
(xj^ and X 2 ) by means of a reduced-order filter. Apply- 
ing the technique discussed in Section A of this chapter, 

t h 

the original model (4^*^order) could .be partitioned into 
two second-order models, but in order to make the compu- 
tational procedures much simpler, the original model was 
partitioned into four scalar models. This partitioning 
approach facilitated the selection of the A^ matrices and 
the computation of their respective pseudo-inverses, A^. 

The solution to the partitioning problem was achiev- 
ed by trial and error. Four new vectors , i = 1,2, 3, 4 
were defined by 



y^(k) = A^x(k) > i = 1,2, 3, 4. 



The four A^ (1x4) matrices selected were 
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0 



0 



0 




A2 = '0 1 0 0 

= 0 0 1 0 

A^ = ,0 0 0 1 

It is important to remember here that all the A^ matrices 
should have the same dimensions, and that the choice of the 
A^*s in this case renders the y^*s scalar quantities. 

The selection of the A^ matrices was made by trial and 
error (in this case obvious and simple) in order to get the 
values of proper A^'s^ that would fulfill the requirement 
established by Eq . (4.3), in this particular case 

4 

^ aIa. = I . 

i. . 11 

i=l 



The pseudo-inverse for each A^ matrix was found using 
Eq. (4.2) 



A. 



1 




A. 



1 



A. 



1 






(4.2) 



1 

Werther (Ref. 9) presents a computation method that 
yields to matrix inversion or pseudo-inversion. 
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from which the four (4 x 1) matrices were found to be^ 



0 




0 



0 



' 0 

I 






The partitioning of the original state vector is then 
summarized by 



yi = : 


1 


0 


0 


0, x(l<) 


= Xj^(k) 


Y2 = 


0 


1 


0 


0|x(k) 


= X2(k) 


ys = 


0 


0 


1 


0"x(k) 


= x^Ck) 


y4 = 


0 


0 


0 


Ij x0<) 


= x^(k) 



with dynamics as specified by Eq. (4.6) 

This partitioning can be interpreted as an apparent 
canonization of the fourth-order model into four first- 
order models, represented each one by one state of the 
original system. 

Having selected the and A^ matrices, four 
coupled scalar filters were implemented in parallel using 



1 

It can be easily verified, that the selection of the A^ 
and A"! matrices agrees with Eqs. (4.2) and (4.3). 
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Eqs. (4.13) and (4.14), but due to the presence of a deter- 
ministic forcing function, Eq . (4.13) was changed slightly 

as follows. 



^.(k/k) = A.$ y At9.(k-l/k-l)+qu(k-l)+G^(k) ^z(k) 



4 

-H , ' At[y.(k/k-l) + nu(k-l)j i = l,2,3,4 

j = l 

( 4 . 22 ) 



to show the effect of u(k) on the prediction step. 

The following dynamics were used for the computation 
of the filter gains and theoretical covariances. 



$1 = AjJa* 



1 = 1 , 2 , 3, 4 



11 0 0 0 



^1 = 






11 



14 



• 22 



33 



41 



44 



1 ^ 

0 I 
0 

0 . 
J 



/l 




1 
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where . ^ and , are scalars. Similarly, 



r 



1 r- 



0 1 0 Oj 1 



4 



14 



' 22 






33 



L ^41 



Uu 





1 

o 




1 




0 




o 

1 



12 / 22 ’ 
I 3 = / 33 ’ 

and 




Also , 



r. = A.r 

1 1 



i = 1,2, 3, 4 



where 
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and similarly, 



72 - 721 

73 = ?31 



and 



"4 



= 

V 41 



Note that the symbols ... ... etc., 

correspond to scalar dynamics of the reduced processes, 

whereas the symbols /n> I 22 ' *** ^11’ ^21* *** > 

correspond to scalar elements of the dynamics of the full 
system. 

The Q matrix was partitioned as follows: 



Q. = A.QAT 



i = 1,2, 3, 4. 



then 



■10 0 0 



Qi = 



'll 



! • 

: ^^41 



• • 



'14 



22 



33 



44 - 

J 



I 

I " 

' 0 

0 ' 



Qi qii , 
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and similarly 



^2 ^22 ’ 
Q 3 = q33 , 



and 



^4 ^44 ’ 

where q 2 ••• q]^]^> ^22 represent scalar quan- 

t it ies . 

The measurement process was specified by 
z(k) = H ^ Aty^(k) + v(k) 

j = l 

with the measurement matrix H = [l 0 0 Oj, and 

v(k) : N(0.0, 0.25). 

For gains and variance calculations for each filter, 
the measurement matrix H(1 x 4) was partitioned as fol- 
lows : 

= HAt i = 1,2, 3, 4 

from which 

h, = ii 0 0 0 1 r 1 ' 

Xu u ; 

0 ; 

0 : 

. 0 J 
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or 






and 




with , h -2 ... being scalars. 

The initialization values for the reduced filters 
were computed using Eqs. (4.15) and (4.16), from which 
the initial guesses for the covariance of estimation 
error P^(0/-1) and expected value of the state vector 
ii(o/-i) were obtained as follows: 



P.(0/-1) = A.P(0/-l)At i = 1,2, 3, 4 



which gives 



[i 0 0 o]rp^^(o/-i) 





r; 



0 




P33(0/-l) 



0 



p^4(0/-l) 0 
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or 



p^(0/-l) = p^^(0/-l) = 0.5 , 
P2(0/-l) = P22(0/-1) = 2.0 , 
P3(0/-l) = P33(0/-1) = 2.0 , 



and 



p^(0/-l) = p^^(0/-l) = 4.0 , 

where p^CO/-!) ..., pj^ 3 ( 0 /-l) ..., are scalars, with 
Pj^(0/-l), p 2 ( 0 /-l), ... etc., representing the initial 
values of covariance for the reduced filters. 

For the initial values of the estimates it is easily 
seen that 

yj^CO/-!) = X3^(0/-l) = 0.0 , 

y2(0/-l) = X2(0/-l) = 0.0 , 

93(07-1) = 93(07-1) = 0.0 , 

and 



9 ^( 07 -!) = 9 ^( 07 -!) = 0.0 

will provide the filters with unbiased initial values. 
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Then four scalar, discrete, sub-opt iraal filters, in 
this example coupled^ only through the measurement process, 
were implemented and run in parallel in accordance with 
Eqs. (4.22), (4.14), and (4.7). At the end of each sam- 
pling interval the outputs from the four filters were com- 
bined by the expression 

4 

x(k/k) = y Aty^(k/k) , 

j=l 

• 

in order to recover the sub-optimal estimates of the 
original vector x(k), from which only the estimation of 
Xj^ and ^2 was of particular interest. 

At the same time, the theoretical sub-optimal covari- 
ance of estimation error ?(k), was computed using recur- 
sive matrix Eq. (4.21), in order to evaluate the selection 
of the set of constant matrices A^. This is possible to 
see when ^(k) is compared with the experimental reduced 
covariance of estimation error. 

2 

Again a Monte Carlo simulation was performed with 
200 replications of the trajectory, each trajectory made 



1 

This is not generally so. The uncoupling here is a con- 
sequence of the choice of the Aj^ matrices. The yj^ state 
equations need not be uncoupled, but the associated gain 
and variance equations are, however, always uncoupled. 

2 

As in Example 1, a computer System 360/67 was used for 
the simulation and solution of this problem. See Appendix 
C for corresponding computer program. 
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of 50 samples, with sampling intervals of 0.05 seconds. As 
a result of this simulation, ensemble values for the mean 
and experimental covariance of estimation error, as a func- 
tion of time, were obtained in and X 2 for the optimal 
and sub-optimal filters as a means of comparing their per- 
formance . 

Numerical results of this simulation are tabulated in 
Tables (lV-2) through (lV-6) and graphical representation 
of the results is shown in Figs. (4.3) through (4.6). 



Table IV-1 



PHI MATRIX 

0.9999 0 

-0.0066 0 

-0.3644 -0 

-12.2713 -10 



0499 0.0012 
9946 0.0480 
2981 0.8871 
1815 -3.9181 



GAMMA MATRIX 

0.0001 

0.0066 

0.3644 

12.2713 

Q MATRIX 

0.0000 0.0000 0.0000 

0.0000 0.0000 0.0000 

0.0000 0.0000 0.0013 

0.0000 0.0008 0.0447 



0 . 0000 
0.0009 
0.0307 
0.3043 



0.0000 

0.0008 

0.0447 

1.5059 



VARIANCE OF FORCING NOISE =0.01 
VARIANCE OF MEASUREMENT NOISE =0.25 



INITIALIZATION - P(0/-1) 



0.5000 

0.0 

0.0 

0.0 



0.0 

2.0000 

0.0 

0.0 



0.0 

0.0 

2.0000 

0.0 



0.0 



0.0 

0.0 

4.0000 



INITIALIZATION - X(0/-1) 

0.0 

0.0 

0.0 

0.0 
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Table IV-2 



GAINS FOR FULL AND REDUCED FILTERS 



K 


G(1 »1 ) 


G(2»l ) 


G(3»l) 


G(4» 1) 


Gl( It 1) 


G2( Itl) 


1 


0.6667 


0.0 


0.0 


0.0 


0.6667 


0.0 


2 


0.4071 


0.2333 


-0.2095 


-7.2845 


0. 4000 


0.0 


3 


0.3097 


0.4075 


-0.6537 


-9.4633 


0.2857 


0.0 


4 


0.2677 


0.5088 


-1.1719 


-9.0573 


0.2222 


0. 0 


5 


0.2477 


0.5314 


-1.6080 


-7.0288 


0.1817 


0.0 


6 


0.2349 


0.4901 


-1.8640 


-4.4015 


0. 1538 


0.0 


7 


0.2227 


0.41Q8 


-1.9266 


-1.9563 


0.1333 


0.0 


8 


0.2087 


0.3175 


-1.8378 


-0.0721 


0.1176 


0.0 


9 


0.1928 


0.2268 


-1.6557 


1. 1851 


0. 1052 


0.0 


10 


0.1758 


0. 1474 


-1.4305 


1.9077 


0.0952 


0.0 


11 


0.1584 


0.0825 


-1.1970 


2.2348 


0.0869 


0.0 


12 


0.1412 


0.0322 


-0.9760 


2.2961 


0. 0799 


0.0 


13 


0.1247 


-0.0051 


-0.7780 


2.1931 


0.0740 


0.0 


14 


0. 1093 


-0.0315 


-0.6072 


1.9984 


0. 0689 


0.0 


15 


0.0950 


-0.0490 


-0.4637 


1.7610 


0.0644 


0.0 


16 


0.0819 


-0.0596 


-0.3457 


1.5120 


0.0605 


0.0 


17 


0.0701 


-0.0647 


-0.2505 


1.2705 


0.0571 


0.0 


18 


0.0596 


-0.0660 


-0.1752 


1.0474 


0.0540 


0.0 


19 


0.0502 


-0.0643 


-0.1167 


0. 8480 


0.0512 


0.0 


20 


0.0420 


-0.0606 


-0.0723 


0.6746 


0.0487 


0.0 


21 


0. 0348 


-0.0556 


-0.0395 


0.5269 


0.0464 


0.0 


22 


0.0286 


-0.0499 


-0.0161 


0.4036 


0. 0444 


0.0 


2.3 


0.0233 


-0.0438 


-0.0001 


0.3026 


0.0425 


0.0 


24 


0. 0188 


-0.0378 


0.0100 


0.2215 


0. 0407 


0.0 


25 


0.0151 


-0.0320 


0.0157 


0. 1576 


0. 0391 


0.0 


26 


0.0120 


-0.0266 


0.0181 


0.1084 


0.0377 


0.0 


27 


0.0095 


-0.0217 


0.0181 


0.0714 


0.0363 


0.0 


28 


0.0075 


-0.0174 


0.0166 


0.0444 


0.0350 


0.0 


29 


0.0059 


-0.0136 


0.0141 


0.0254 


0.0338 


0.0 


30 


0.0047 


-0.0105 


0.0111 


0.0127 


0.0327 


0.0 


31 


0.0037 


-0.0078 


0.0080 


0.0046 


0.0317 


0.0 


32 


0. 0031 


-0.0057 


0.0050 


-0.0000 


0.0307 


0.0 . 


33 


0.0026 


-0.0041 


0.0023 


-0.0022 


0. 0298 


0.0 


34 


0.0022 


-0.0028 


-0.0001 


-0.0028 


0.0289 


0.0 


35 


0.0020 


-0.0018 


-0.0022 


-0.0024 


0.0281 


0.0 


36 


0.0018 


-0.0011 


-0.0038 


-0.0015 


0.0273 


0.0 


37 


0. 0018 


-0.0006 


-0.0050 


-0.0004 


0.0266 


0.0 


38 


0.0017 


-0.0003 


-0.0060 


0.0007 


0.0259 


0.0 


39 


0.0017 


-0.0001 


-0.0066 


0.0017 


0.0252 


0.0 


40 


0.0017 


0.0000 


-0.0070 


0.0024 


0.0246 


0.0 


41 


0.0017 


0.0001 


-0.0073 


0.0030 


0. 0240 


0.0 


42 


0.0017 


0.0001 


-0.0074 


0.0033 


0.0234 


0.0 


43 


0.0017 


0.0001 


-0.0074 


0.0034 


0.0229 


0.0 


44 


0.0017 


0.0001 


-0.0073 


0.0033 


0.0224 


0.0 


45 


0.0017 


0.0000 


-0.0072 


0.0031 


0.0219 


0.0 


46 


0.0017 


0.0000 


-0.0071 


0.0028 


0.0214 


0.0 


47 


0.0017 


-0.0000 


-0.0070 


0.0025 


0.0210 


0.0 


48 


0.0017 


-0. 0000 


-0.0068 


0.0022 


0.0205 


0.0 


49 


0.0017 


-0.0001 


-0.0067 


0.0018 


0.0201 


0.0 


50 


0. 0017 


-0.0001 


-0.0066 


0.0015 


0.0197 


0.0 

i 




G (4x1) 


represents the matrix of gains for the 


optimal 



filter. 



G1 and G2 represent the gains used in the sub-optimal 
filters 1 and 2. 
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Table IV-3 



THEORETICAL 
FULL PLANT 



VARIANCE OF ESTIMATION ERROR 



K 


PKI 1, 1) 


PK(2f2l 


PK ( 3 , 3 ) 


PK(4,4) 


1 


0.1667 


2.0000 


2.0000 


4.0000 


2 


0.1018 


1. 9603 


1.7605 


242.6279 


3 


0.0774 


1.8002 


2.3189 


346.3362 


4 


0.0669 


1.5182 


3.9792 


261.1479 


5 


0.0619 


1.1771 


5.7036 


142.2309 


6 


0.0587 


0.8478 


6.6569 


60.6997 


7 


0.0557 


0.5745 


6.6825 


21.320Q 


8 


0.0522 


0.3709 


6.0419 


8.2223 


9 


0.0482 


0.2307 


5.0782 


6.9855 


10 


0.0439 


0. 1400 


4.0514 


9.3412 


11 


0.0396 


0.0847 


3.1097 


11.6978 


12 


0.0353 


0.0529 


2.3160 


12.9672 


13 


0.0312 


0.0361 


1.6824 


13. 1242 


14 


0.0273 


0.0282 


1.1960 


12.4778 


15 


0.0237 


0. 0255 


0.8337 


11.3650 


16 


0.0205 


0.0254 


0.5710 


10.0516 


17 


0.0175 


0.0264 


0.3855 


8.7171 


18 


0. 0149 


0.0275 


0.2581 


7,4696 


19 


0.0126 


0.0283 


0.1736 


6, 3652 


20 


0.0105 


0. 0286 


0.1199 


5.4260 


21 


0.0087 


0.0282 


0.0877 


4,6526 


22 


0.0072 


0,0273 


0.0701 


4.0336 


23 


0. 0058 


0.0259 


0.0619 


3.5512 


24 


0.0047 


0.0241 


0.0595 


3. 1854 


25 


0. 0038 


0.0221 


0.0601 


2.9159 


26 


0. 0030 


0,0199 


0.0619 


2.7237 


27 


0.0024 


0.0177 


0.0639 


2. 5918 


28 


0.0019 


0, 0155 


0.0654 


2.5057 


29 


0.0015 


0.0133 


0. 0659 


2.4529 


30 


0.0012 


0.0114 


0.0654 


2.4236 


31 


0. 0009 


0.0096 


0.0640 


2.4099 


32 


0.0008 


0.0081 


0.0618 


2.4059 


33 


0.0006 


0.0067 


C.0590 


2.4071 


34 


0. 0006 


0.0056 


0.0557 


2.4105 


35 


0.0005 


0.0046 


0.0523 


2.4141 


36 


0.0005 


0.0039 


0.0489 


2.4167 


37 


0.0004 


0.0033 


0,0455 


2.4177 


38 


0.0004 


0.0028 


0.0424 


2.4170 


39 


0. 0004 


0.0024 


0.0396 


2.4146 


40 


0.0004 


0.0022 


0.0371 


2.4109 


41 


0.0004 


0.0020 


0.0349 


2.4062 


42 


0.0004 


0.0018 


0.0330 


2.4008 


43 


0.0004 


0,0017 


0.0315 


2. 3951 


44 


0.0004 


0.0017 


0.0303 


2.3894 


45 


0.0004 


0.0017 


0.0293 


2.3839 


46 


0.0004 


0.0016 


0.0285 


2.3787 


47 


0.0004 


0.0016 


0.0279 


2.3741 


48 


0.0004 


0.0016 


0.0275 


2.3700 


49 


0.0004 


0.0016 


0.0272 


2.3664 


50 


0.0004 


0.0016 


0.0270 


2.3635 



PK represents 
covariance 



the diagonal elements of the theoretical 
matrix for the optimal filter. 
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Table IV-4 



EXPERIMENTAL VARIANCE OF ESTIMATION ERROR 
FULL PLANT 



K 


VKI 1, 1) 


VKI 2 


>2) 


VK(3»3) 


VKI4,4I 


1 


0.1317 


1.9066 


1.9733 


4. 0880 


2 


0.0910 


1.8538 


1.7402 


225.0217 


3 


0.0739 


1.7340 


2.2403 


326.4661 


4 


0.0623 


1.4891 


3.7969 


257. 7678 


5 


0.0645 


1.1562 


5.6106 


138.7382 


6 


0.0518 


0.8086 


5.9918 


60.5067 


7 


0.0510 


0.5760 


6.5021 


23.3325 


8 


0.0485 


0.3807 


6.0129 


8.6785 


9 


0.0465 


0.2391 


5.1270 


7.2472 


10 


0.0459 


0.1500 


4.3535 


9.5983 


11 


0.0424 


0.0927 


3.5224 


12.3814 


12 


0.0394 


0.0557 


2.7024 


14.0400 


13 


0.0367 


0.0360 


1.9633 


13.9882 


14 


0.0295 


0.0276 


1.3043 


13.4596 


15 


0.0264 


0.0240 


0.9568 


11.7134 


16 


0.0225 


0.0239 


0.6456 


12.1689 


17 


0.0187 


0.0256 


0.4017 


9.6542 


18 


0.0152 


o.o: 




0.2554 


7.2280 


19 


0.0130 


o.oi 


>71 


0.1698 


5.9514 


20 


0.0111 


0.0269 


0.1194 


5.5124 


21 


0.0093 


0.0275 


0.0835 


5.0569 


22 


0.0073 


0.0253 


0.0671 


4.1066 


23 


0.0060 


0.0241 


0.0626 


3.5153 


24 


0.0047 


0.0218 


0.0616 


2.8292 


25 


0.0039 


0.0203 


0.0586 


3.1777 


26 


0.0031 


0.0186 


0.0564 


2.4777 


27 


0.0025 


0.0] 


L66 


0.0546 


2.7436 


28 


0.0020 


0.0 


150 


0.0546 


2.7655 


29 


0.0016 


0.0] 


[31 


0.0550 


2.6672 


30 


0.0012 


0.01 


L15 


0.0581 


2.0752 


31 


0.0010 


0.0099 


0.0586 


2.4300 


32 


0.0008 


0.0084 


0.0572 


2. 7947 


33 


0.0007 


0.0072 


0.0546 


2.5196 


34 


0.0006 


0.0060 


0.0495 


2.1787 


35 


0.0005 


0.0051 


0.0437 


2.5699 


36 


0.0004 


0.0043 


0.0422 


2.2536 


37 


0.0004 


0.0035 


0.0438 


2.5338 


38 


0.0004 


0.0029 


0.0449 


2.3031 


39 


0.0004 


0.0025 


0.0427 


2.0479 


40 


0.0004 


0.0021 


0.0373 


2.4109 


41 


0.0004 


0.0019 


0.0308 


2.2848 


42 


0.0004 


0.0018 


0.0261 


2.5309 


43 


0.0004 


0.0016 


0.0242 


2.2161 


44 


0.0004 


0.0015 


0.0261 


2.4027 


45 


0.0004 


0.0014 


0.0299 


1.9605 


46 


0.0004 


0.0014 


0.0316 


1.9165 


47 


0.0004 


0.0014 


0.0306 


2.2086 


48 


0.0004 


0.0014 


0.0297 


2.3207 


49 


0.0004 


0.0015 


0.0281 


2.1990 


50 


0.0004 


0.0016 


0.0229 


2.8940 



VK represents the diagonal elements of the ex- 
perimental covariance matrix for the optimal 
process . 
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Table IV-5 



SUBOPTIMAL VARIANCE OF ESTIMATION ERROR 



K 


SR(1,1) 


SR<2,2) 


VR<1,U 


VR(2t2l 


1 


0.1667 


2.0000 


0.1308 


1.9041 


2 


0.1018 


1.9833 


0.0899 


1.8888 


3 


0.0777 


1. 8799 


0.0725 


1.7924 


4 


0.0682 


1.6698 


0.0615 


1.5970 


5 


0.0653 


1.3827 


0.0627 


1.3288 


6 


0.0653 


1.0666 


0.0578 


1.0320 


7 


0.0661 


0.7655 


0.0603 


0.7471 


8 


0.0662 


0.5098 


0. 0610 


0. 5028 


9 


0.0650 


0.3136 


0.0603 


0.3134 


10 


0.0622 


0.1780 


0.0598 


0.1806 


11 


0.0579 


0.0954 


0.0553 


0.0990 


12 


0.0525 


0.0538 


0.0532 


0.0549 


13 


0.0464 


0. 0406 


0.0490 


0.0393 


14 


0.0399 


0.0443 


0.0406 


0.0407 


15 


0.0336 


0.0558 


0.0351 


0.0506 


16 


0.0277 


0.0687 


0. 0294 


0.0627 


17 


0.0225 


0.0791 


0.0241 


0.0729 


18 


0.0180 


0.0852 


0.0187 


0.0793 


19 


0.0143 


0.0864 


0.0149 


0.0813 


20 


0.0114 


0.0834 


0.0120 


0.0794 


21 


0.0091 


0. 0770 


0.0096 


0.0743 


22 


0.0075 


0.0685 


0.0072 


0.0669 


23 


0.0063 


0.0589 


0.0061 


0.0582 


24 


0.0055 


0.0492 


0.0048 


0.0493 


25 


0.0049 


0.0401 


0.0046 


0.0406 


26 


0.0046 


0.0320 


0.0041 


0.0327 


27 


0.0044 


0.0251 


0.0038 


0.0258 


28 


0.0042 


0.0196 


0.0040 


0.0202 


29 


0.0041 


0.0153 


0.0039 


0.0157 


30 


0.0039 


0.0121 


0.0038 


0.0123 


31 


0.0038 


0.0099 


0.0036 


0.0099 


32 


0.0037 


0.0084 


0.0035 


0.0083 


33 


0.0035 


0.0074 


0.0034 


0.0074 


34 


0.0033 


0.0068 


0.0030 


0.0069 


35 


0.0032 


0.0065 


0.0028 


0.0066 


36 


0.0030 


0.0063 


0.0026 


0.0063 


37 


0.0028 


0.0062 


0.0025 


0.0061 


38 


0.0026 


0.0061 


0.0025 


0.0058 


39 


0.0025 


0.0060 


0.0026 


0.0056 


40 


0.0023 


0.0058 


0.0026 


0.0054 


41 


0.0022 


0.0056 


0.0025 


0.0052 


42 


0.0021 


0.0054 


0.0023 


0.0050 


43 


0.0020 


0.0052 


0.0023 


0.0047 


44 


0.0019 


0.0050 


0.0022 


0.0045 


45 


0.0018 


0. 0047 


0.0020 


0.0042 


46 


0.0017 


0.0045 


0.0020 


0.0040 


47 


0.0016 


0.0043 


0.0018 


0.0039 


48 


0.0016 


0.0041 


0.0017 


0.0038 


49 


0.0015 


0.0039 


0.0017 


0.0038 


50 


0.0015 


0.0037 


0.0017 


0.0038 



SR represents two of the diagonal elements of the sub 
optimal theoretical covariance of estimation error 
corresponding to the sub-optimal estimation of x, 
and X 2 > 

VR follows exactly the same pattern as SR, but it re- 
fers to the experimental covariance of estimation 
error. 
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Table IV-6 



EXPERIMENTAL MEANS OF ESTIMATION ERROR 



K 


MEANEl 


MEANE2 


MEANE3 


MEANE4 


MEANEIR 


MEANE2R 


1 


0.0312 


0.0500 


-0.0350 


-0. 1543 


0.0312 


0.0500 


2 


0.0263 


0.0436 


-0.0543 


-0.5355 


0.0264 


0.0478 


3 


0.0282 


0.0398 


-0. 0893 


-0.7945 


0.0285 


0.0439 


4 


0.0330 


0.0398 


-0.1411 


-0. 8346 


0.0328 


0.0379 


5 


0.0362 


0.0349 


-0.1811 


-0.4616 


0. 0356 


0.0301 


6 


0.0313 


0.0120 


-0. 1482 


-0. 1316 


0. 0328 


0.0213 


7 


0.0289 


-0.0006 


-0.1229 


0. 1844 


0.0317 


0.0120 


8 


0.0041 


-0.0439 


0. 1008 


0.1273 


0.0178 


0.0030 


9 


-0.0014 


-0.0429 


0.1368 


0.0809 


0. 0142 


-0.0054 


10 


0.0055 


-0.0284 


0.0687 


0. 1810 


0. 0170 


-0.0124 


11 


0.0150 


-0.0193 


-0.0130 


0.0214 


0.0211 


-0.0176 


12 


0.0031 


-0.0223 


0. 0684 


-0.0040 


0. 0134 


-0.0213 


13 


-0.0017 


-0.0188 


0.0900 


-0. 1130 


0. 0093 


-0.0238 


14 


-0. 0060 


-0.0134 


0.1014 


-0.2676 


0. 0052 


-0.0251 


15 


-0.0075 


-0.0083 


0.0901 


-0.3582 


0. 0026 


-0. 0254 


16 


-0.0046 


-0.0065 


0.0577 


-0.3345 


0.0031 


-0.0250 


17 


-0.0016 


-0.0069 


0.0383 


0.0395 


0.0042 


-0.0239 


18 


0.0031 


-0.0106 


0.0206 


-0.0433 


0. 0073 


-0.0220 


19 


0.0066 


-0.0148 


0.0064 


-0.0677 


0.0101 


-0.0198 


20 


0.0054 


-0.0139 


0.0036 


-0.0764 


0. 0085 


-0.0179 


21 


0.0050 


-0.0140 


0.0084 


0.2286 


0.0078 


-0.0164 


22 


0.0024 


-0. 0100 


0.0190 


0.1291 


0.0040 


-0.0148 


23 


0.0030 


-0.0111 


0.0238 


0.0811 


0.0053 


-0.0131 


24 


0.0025 


-0. 0099 


0.0236 


-0. 0649 


0.0046 


-0.0115 


25 


0.0018 


-0.0082 


0.0217 


-0.0134 


0. 0033 


-0. 0102 


26 


0.0015 


-0.0074 


0.0204 


-0.0375 


0.0030 


-0. 0090 


27 


0.0013 


-0.0068 


0.0202 


0.0076 


0.0032 


-0.0080 


28 


0.0011 


-0.0059 


0.0262 


0. 1873 


0. 0034 


-0. 0069 


29 


0.0006 


-0.0038 


0.0332 


0. 1184 


0.0015 


-0.0055 


30 


0. 0006 


-0.0025 


0.0299 


-0.2134 


0.0022 


-0.0040 


31 


0.0007 


-0.0016 


0.0225 


-0. 1106 


0. 0036 


-0.0028 


32 


0.0007 


-0.0008 


0.0181 


-0.0725 


0.0043 


-0.0019 


33 


0.0008 


-0.0001 


0.0093 


-0.2493 


0. 0044 


-0.0013 


34 


0.0007 


0.0003 


0.0026 


-0.0454 


0.0033 


-0.0013 


35 


0. 0006 


0.0005 


0.0036 


0.0645 


0.0020 


-0.0014 


36 


0.0007 


0.0006 


0.0030 


-0.0675 


0. 0020 


-0.0015 


37 


0.0007 


0.0007 


0.0003 


-0.0409 


0.0020 


-0.0016 


38 


0. 0007 


0.0006 


-0.0055 


-0. 1695 


0.0014 


-0.0019 


39 


0.0007 


0.0002 


-0.0080 


0.0391 


0.0012 


-0.0024 


40 


0.0007 


-0.0001 


-0.0045 


0.0863 


0.0005 


-0.0027 


41 


0.0008 


-0.0002 


0.0006 


0. 1190 


0.0015 


-0.0027 


42 


0.0007 


-0.0001 


0.0038 


0.0133 


0.0005 


-0.0023 


43 


0.0008 


0.0002 


0.0048 


0.0356 


0.0015 


-0.0019 


44 


0.0008 


0.0004 


0.0035 


-0. 0646 


0. 0019 


-0.0014 


45 


0.0008 


0.0005 


0.0012 


-0.0425 


0.0009 


-0.0010 


46 


0.0008 


0.0006 


0.0043 


0. 1323 


0. 0006 


-0.0007 


47 


0.0007 


0.0009 


0.0081 


0.0145 


-0.0009 


-0.0002 


48 


0.0008 


0.0013 


0.0074 


-0.0337 


-0.0008 


0.0003 


49 


0.0009 


0.0017 


0.0071 


0.0159 


-0.0005 


0.0009 


50 


0.0010 


0.002C 


0.0057 


-0.0529 


0. 0001 


0.0015 



MEANEl through MEANE4 represent the mean of estimation 
error in the complete state vector by means of optimal 
filtering . 

MEANEIR and MEANE2R represent the sub-optimal mean of 
estimation error in states x-i and Xo by means of reduced 
filters . 
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Variance of estimation error 




Fig. 4.3 - Graphical plot of the variance of estima- 
tion error in xi vs. time. Curves a and b show the 
experimental and theoretical variance of estimation 
error for the optimal process. Curves c and d show 
the sub-optimal experimental and theoretical co- 
variance of estimation error. 
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Variance of estimation error 




Fig. 4.4 - Graphical plot of the variance of estima- 
tion error in X2 vs. time. Curves a and b show the 
experimental and theoretical variance of estimation 
error for the optimal process. Curves c and d show 
the sub-optimal case. 
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Mean of estimation error 




Fig. 4.5 - Graphical plot of the mean of estimation 
error in vs. time. Curve a shows the mean of 
estimation error for the optimal estimation process. 
Curve b shows that for the sub-optimal case. 
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Mean of estimation error 



0.08 




Time in seconds 

Fig. 4.6 - Graphical plot of the means of estimation 
error in X2 vs. time. Curve a shows the mean of 
estimation error for the optimal estimation process. 
Curve b shows that for the sub-optimal case. 
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V. IMPLEMENTATION OF A REDUCED DISCRETE 



KALMAN FILTER BY USING THE DOMINANT 
ROOTS OF THE ORIGINAL PROCESS 



The general idea of producing a reduced sub-optimal 
filter (Chapter III) or several sub-optimal filters of re- 
duced order (Chapter IV), is presented here using a dif- 
ferent approach. 

A. THE DOMINANT-ROOTS APPROACH 

It is assumed as in the previous chapters, that the 
full system dynamics are represented by 

x(k+I) = r j^u(k) + w(k)j , (2.9) 

with measurement 

^(k) = Hx(k) + v(k) , (2.3) 

and with the random processes involved, w(k) and v(k), 
fully described by Eqs. (2. II) through (2.16). 

It is clear that the order of the system is n and 
that the goal of the problem is to estimate only r states 
(r<n) of the process; therefore it will be necessary to 
produce a reduced model of the plant in order to use its 
dynamics in the implementation of the reduced filter. 
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It is common practice when handling classical auto- 
matic control problems to assume, or to require, that the 
high-order system be dominated by one pair of complex con- 
jugate roots. This will normally permit the use of a se- 
cond-order characteristic equation, and will simplify 
analysis and synthesis because the second-order system has 

been completely solved and all kind of curves, etc., are 

t ll 

available. Applying the same idea to the n^ -order, linear 
plant under consideration, but without the constraint of a 
pure second-order approximation, the dominant -root s repre- 
sentation will be determined mainly by the number r of 
states to be estimated and by the amount of performance de- 
gradation introduced when choosing the reduced model. 

It is noted that the number r of states to be estimat- 
ed will specify the minimum number of roots which can be 
selected as dominant, but an upper limit in the selection 
of the dominant roots would be determined by the designer's 
judgement, based upon the amount of degradation introduced 
by the reduced model. 

Due to the qualitative nature of the root -dominance 
criterion, it is very difficult to state a general proce- 
dure for the selection of the reduced-order model. There- 
fore a specific example will be used to illustrate the 
method . 
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Consider the fourth-order linear system with continu- 
ous dynamics used in Examples 1 and 2 and described by 



A(t) = 



r 0 



0 

0 

-400 



1 

0 

0 

-320 



0 0 

1 0 

0 i! 

•118 -19i 



0 ^ |'u(t)+w(t) 



x(t) + 



0 

0 

;400 



J 



with I.C 



x(0) = 0 , 



and measurement 



(5.1) 



;(t) = 



10 0 0 



x(t) + v(t) 



The task of the problem is to estimate only two states 
of the original process (xj^ and X 2 ). 

When solving Example 1 (Chapter 3), the whole process 
(Eq. 5.1) was discretized and then the dynamics of the dis- 
crete system were partitioned according to the number of 
states to be estimated. Here, instead of discretizing the 
whole system, the order of the continuous process will be 
lowered by proper choice of the dominant roots, and then 
the reduced continuous system will be discretized. 
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The system described by Eq . (5.1) may be rewritten as 



Xj^(t) = X2(t) , 

X2(t) = x^Ct) , 
X 2 (t) = x^(t) , 

and 



(5.2a) 

(5.2b) 

(5.2c) 



^4(t) = -400x^ ( t ) -320x2 ( t )-118 x 2( t ) -19x^( t )+400 u(t)+w(t) . 

(5. 2d) 

For convenience a new variable .^(t) = x^(t) is defined; 
then Eqs. (5.2a,b,c) become 



“(t) = X2(t) , 
t7(t) = x^(t) , 
t7(t) = x^(t) , 

and by direct substitution of these three equations into 
Eq . (5. 2d) the latter becomes 



C^(t) + 19 9 (t) + 118 (9 (t) + 320e(t) +400-;'(t) =400[u(t)+w(t)j , 



(5.3) 



where 



(t) = 






dt 
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The dominant-root analysis of Eq. (5.3), which now fully 
represents the fourth-order system under consideration, 
will be made with the help of a pole-zero representation. 
This will require (as an easy approach) transformation of 
Eq. (5.3) into the s-domain, accomplished here by use of 
the Laplace transform. Then, considering all initial con- 
ditions equal to zero, Eq . (5.3) transforms into 

®(s)s^+19fl)(s)s^+1180(s)s^+32O®(s)s+4OO0(s)=4OO [u(s)+W(s)] . 



which after proper factorization gives 



®(s) = 



400.U(s)+W(s) 



] 



(s+10)(s+5)(s+2-j2)(s+2+j2) 



or 



®(s) ^ 400 

U(s)+W(s) (s+10)(s+5)(s+2-j2)(s+2+j2) 

(5.4) 

Eq. (5.4), may be interpreted as the transfer function of 
the fourth-order system relating (t$(s) with the forcing 
functions (inputs). 

A pole-zero representation of Eq . (5.4) (Fig. 5. lb), 
shows a pair of complex conjugate roots located at 
s = -2+j2, and two real roots located at s = -5 and s = -10 
respectively. If the dominance criterion is based upon the 
size of the real portion of the roots under consideration. 
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Fig. 5.1 - a). Block diagram representation of the pro- 
cess described by Eq. (5.4). 

b). Pole-zero diagram for Eq. (5.4) indicating four roots 
s = -5, s *= “10, s “ -2ij2 respectively, and with 400 
considered to be the gain of the system. 
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then it is clear for this particular example, that the less 
dominant root will be at s = -10 (this will give rise to a 
term of the form £ and the most dominant roots will 

be those at s = -2+j2, with the root at s = -5 producing an 
intermediate effect. 

As stated before, the purpose of this problem is to 
estimate only two states of the original system (xj^ and X 2 ). 
It therefore seems a good approach to select the roots at 
s = -2ij2 in order to produce a reduced second-order con- 
tinuous system, whose dynamics will be used in the imple- 
mentation of the reduced filter. It is important to note 
that the gain of the reduced system should be corrected in 
order to compensate for the disappearance of two of the 
roots of the original system. In this case the gain of the 
reduced system was made equal to the product of the magni- 
tudes of the dominant roots, that is 



gain reduction 



-2-j2| j-2+j2 



8 . 



See Fig. (5.2) 

Another important consideration should be based upon 
the fact that the performance degradation introduced, and 
the number m of inputs (forcing functions) to the original 
system, will strongly determine the order of the reduced 
model when using dominant roots. 
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W(s) 




a) . 




Fig. 5.2 - a). Block diagram representation of the 
reduced model. 

b) . Pole-zero diagram showing only the dominant 
roots. Gain of the system has been adjusted to 8. 
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After the selection of the pair of complex conjugate 
roots at s = -2jlj2, the reduced model may be expressed as 



^ 8 ^ 8 

U(s)+W(s) (s+2-j2) (s+2+j2) s2+4s+8 

or transferring back to the differential equation form, 

+ 4^(t) + S6^(t) - 8[u(t)+w(t)] . 

Recalling that A(t) = x^(t), the second-order reduced 
model becomes 



Xj^(t) = X2(t) 

X 2 (t) = -8xj^(t) - 4 x 2 ( 1 ) + 8 j^u(t)+w(t)j , 
or 



X (t) 
r ' 



0 


-1 




0 


1 




Xj.(t) + 




[-8 


-4 




8 



j^u(t)+w(t) 1 



J 



which represents the reduced model from which the reduced 
filter will be implemented (see Fig. 5.2). 

No special derivation of the reduced-filter equations 
is necessary because the form of the equations will be ex- 
actly the same as those for the optimal filtering, and 
only the order will change. 
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The discrete reduced model is obtained from Eq . (5.5) 



using 




A T 



r 



(5,6) 



and 



r 



r 




-A r 

€ ^ dTB , 

r ’ 



(5.7) 



with A^(r X r) and B^(r x ro) representing the reduced sys- 
tem and reduced distribution matrices respectively for the 
continuous case, and f^(r x r) and r^(i^ x m) being the 
state transition matrix and distribution matrix for the 
reduced discrete case. Then the discretized reduced pro- 
cess can be described as 



Xj.(l<+1) = r^j^u(k)+w(k)j 

with measurement 



z(k) = Hx(k) + v(k) . 



As in the partitioned method presented in Chapter 3, 
all the measured outputs are considered to be included in 
the reduced model, and furthermore the H matrix is assumed 
to be partitioned as 



H 



- [HiI H2) , 



with H 2 = 0, and Hj^(p x r) 
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Considering the problem of initialization for the dis- 
crete reduced covariance equations, it is solved exactly as 
in Chapter 3, by partitioning the P(0/-1) matrix from the 
filter, as 



P(0/-1) 



fp . (0/-1) 1 

1 




P2i(0/-1) 1 





and choosing 

P^(0/-1) = P^^(0/-1) 

with Pj^j^(0/-1) (r X r) and P^(0/-1) (r x r) , the latter 
representing the initial value of covariance for the re- 
duced filter. 

For the reduced filter initial states, 
x^(0/-l) = 



where 



x(0/-l) 



Xj^(0/-1) 

X2(0/-l) 
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B. EXAMPLE 3 

The same example used in the previous methods is em- 
ployed here to illustrate the method of dominant roots. 

The same initialization values, same statistics for 
the random sequences, and same measurements as those used 
in Example 1, were used in this example for both the full 
and reduced systems. The corresponding values for the re- 
duced filter are listed here for convenience 

0.5 0 

0 2 . 0 _ 

H = U Oj , , 

and 



0 ! 

Xj,(o/-i) = : , Pj.(o/-i) = 

' 0 ' 



v(k) : N(0.0, 0.25) , 
w(k) : N(0.0, 0.01) . 



A significant difference between the two examples is the 
determination of the reduced transition and distribution 
matrices which in this case were computed as indicated by 
Eqs. (5.6) and (5.7), with 



” 0 


l' 




|0^ 


oo 

L, 


1 

1 


, B 

’ r 


00 
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Due to the similarity between the partitioning method and 
the dominant-roots approach for this particular problem, 
the same computer program was used for the Monte Carlo simula- 
tion and solution of Examples 1 and 3, with the only change 
introduced regarding the computation of the dynamics for 
the reduced model. 

Table (V-l) shows the initialization values for the 
reduced model, as well as the values for the and ma- 
trices. Tables (V-2) through (V-5) present numerical re- 
sults of the simulation, and a graphical representation of 
the simulation is shown in Figs. (5.3) through (5.6). 

Fig. (5.7) shows a complete block diagram of the re- 
duced filtering process for this particular example. 
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Table V-1 



PHI MATRIX 



0.9999 0.0499 0.0012 

-0.0066 0.9946 0.0480 

-0.3644 -0.2981 0.8871 

-12.2713 -10.1815 -3.9181 



0.0000 

0.0009 

0.0307 

0.3043 



GAMMA MATRIX 

0.0001 

0.0066 

0.3644 

12.2713 

REDUCED PHI MATRIX 

0.9906 0.0452 

-0.3613 0.8100 

REDUCED GAMMA MATRIX 

0.0093 

0.3613 

Q MATRIX (REDUCED) 

0.0000 0.0000 

0.0000 0.0000 

VARIANCE OF FORCING NOISE =0.01 

VARIANCE OF MEASUREMENT NOISE =0.25 



REDUCED P(0/-1) 

0.5000 0.0 

0.0 2.0000 

INITIALIZATION - X(0/-1) 



0.0 

0.0 

0.0 

0.0 
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Table V-2 



GAINS FOR FULL AND REDUCED FILTERS 



K 


G ( Itl ) 


G ( 2 , 1 ) 


GI 3 , 1 ) 


G ( 4t 1 ) 


GR ( 1 , 1 ) 


GRI2 , 1 ) 


1 


0.6667 


0.0 


0.0 


0.0 


0 . 6667 


0.0 


2 


0.4071 


0.2333 


- 0.2095 


- 7.2845 


0 . 4014 


0.0323 


3 


0.3097 


0.4075 


- 0.6537 


- 9.4633 


0. 2896 


0.0546 


4 


0.2677 


0.5088 


- 1.1719 


- 9.0573 


0 . 2286 


0.0527 


5 


0.2477 


0.5314 


- 1.6080 


- 7.0288 


0. 1895 


0.0360 


6 


0.2349 


0 . 4901 


- 1.8640 


- 4.4015 


0. 1613 


0.0134 


7 


0.2227 


0.4108 


- 1.9266 


- 1.9563 


0. 1390 


- 0. 0091 


8 


0.2087 


0 . 3175 


- 1.8378 


- 0.0721 


0. 1205 


- 0.0286 




0. 1928 


0.2268 


- 1.6557 


1 . 1851 


0. 1045 


- 0.0436 


10 


0. 1758 


0.1474 


- 1.4305 


1.9077 


0.0903 


- 0.0541 


11 


0 . 15 84 


0.0825 


- 1.1970 


2. 2348 


0.0778 


- 0.0605 


12 


0.1412 


0.0322 


- 0.9760 


2.2961 


0 . 0666 


- 0.0633 


13 


0. 1247 


- 0.0051 


- 0.7780 


2.1931 


0.0566 


- 0.0633 


14 


0. 1093 


- 0.0315 


- 0.6072 


1.9984 


0.0478 


- 0.0612 


15 


0.0950 


- 0.0490 


- 0.4637 


1.7610 


0.0401 


- 0.0576 


16 


0.0819 


- 0.0596 


- 0.3457 


1.5120 


0.0333 


- 0.0529 


17 


0.0701 


- 0.0647 


- 0.2505 


1.2705 


0.0274 


- 0.0477 


18 


0.0596 


- 0.0660 


- 0.1752 


1.0474 


0.0224 


- 0. 0422 


19 


0 . 05 02 


- 0.0643 


- 0.1167 


0.8480 


0.0181 


- 0.0367 


20 


0.0420 


- 0.0606 


- 0.0723 


0.6746 


0. 0144 


- 0.0314 


21 


0. 0348 


- 0.0556 


- 0.0395 


0.5269 


0.0114 


- 0.0264 


22 


0. 0286 


- 0.0499 


- 0.0161 


0.4036 


0.0089 


- 0.0219 


23 


0.0233 


- 0.0438 


- 0.0001 


0.3026 


0. 0069 


- 0.0179 


24 


0.0188 


- 0.0378 


0.0100 


0.2215 


0.0052 


- 0.0143 


25 


0.0151 


- 0.0320 


0. 0157 


0. 1576 


0.0039 


- 0.0113 


26 


0.0120 


- 0.0266 


0.0181 


0. 1084 


0.0029 


- 0. 0083 


27 


0 . 0095 


- 0.0217 


0.0181 


0.0714 


0.0021 


- 0.0067 


28 


0.0075 


- 0.0174 


0.0166 


0.0444 


0.0016 


- 0.0050 


29 


0.0059 


- 0.0136 


0.0141 


0.0254 


0.0011 


- 0.0037 


30 


0.0047 


- 0.0105 


0.0111 


0 . 0127 


0.0008 


- 0.0026 


31 


0.0037 


- 0.0078 


0.0080 


0.0046 


0. 0006 


- 0 . 0018 


32 


0. 0031 


- 0.0057 


0.0050 


- 0.0000 


0.0004 


- 0.0013 


33 


0.0026 


- 0.0041 


0.0023 


- 0. 0022 


0.0003 


- 0.0008 


34 


0.0022 


- 0.0028 


- 0.0001 


- 0.0028 


0. 0003 


- 0.0006 


35 


0. 0020 


- 0.0018 


- 0.0022 


- 0.0024 


0.0002 


- 0.0004 


36 


0.0018 


- 0.0011 


- 0.0038 


- 0.0015 


0. 0002 


- 0. 0002 


37 


0.0018 


- 0.0006 


- 0.0050 


- 0.0004 


0.0002 


- 0.0002 


38 


0.0017 


- 0.0003 


- 0.0060 


0.0007 


0. 0001 


- 0.0001 


39 


0.0017 


- 0.0001 


- 0.0066 


0.0017 


0.0001 


- 0. 0001 


40 


0.0017 


0 . 0000 


- 0.0070 


0.0024 


0. 0001 


- 0.0001 


41 


0.0017 


0 . 0001 


- 0.0073 


0.0030 


0. 0001 


- 0.0001 


42 


0.0017 


0.0001 


- 0.0074 


0.0033 


0. 0001 


- 0.0001 


43 


0.0017 


0.0001 


- 0.0074 


0.0034 


0.0001 


- 0.0001 


44 


0.0017 


0.0001 


- 0.0073 


0.0033 


0. 0001 


- 0.0001 


45 


0.0017 


0.0000 


- 0.0072 


0.0031 


0. 0001 


- 0.0001 


46 


0.0017 


0.0000 


- 0.0071 


0.0028 


0.0001 


- 0.0001 


47 


0.0017 


- 0.0000 


- 0.0070 


0.0025 


0 . 0001 


- 0. 0001 


48 


0.0017 


- 0.0000 


- 0.0068 


0.0022 


0.0001 


- 0.0001 


49 


0.0017 


- 0.0001 


- 0.0067 


0.0018 


0.0000 


- 0.0001 


50 


0.0017 


- 0.0001 


- 0.0066 


0.0015 


0. 0000 


- 0. 0001 




G (4x1) 


represents 

filter. 


the matrix 


of gains 


for the 


opt imal 




GR(2xl) 


represents 

filter. 


the matrix 


of gains 


for the 


reduced 
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Table V-3 



THEORETICAL VARIANCE OF ESTIMATION ERROR 



K 


PK ( l . l ) 


PKI2f2 ) 


PK ( 3,3 ) 


PKI4 , 4 ) 


PRI 1 . 1 ) 


P R ( 2 , 2 ) 


1 


0.1667 


2.0000 


2.0000 


4.0000 


0. 1667 


2. 0000 


2 


0 . 1018 


1. 9603 


1.7605 


242.6279 


0. 1004 


1.3335 


3 


0.0774 


1.8002 


2 . 3189 


346. 3362 


0. 0724 


0.8822 


4 


0.0669 


1.5182 


3.9792 


261.1479 


0.0571 


0. 5793 


5 


0.0619 


1.1771 


5.7036 


142.2309 


0.0474 


0.3794 


6 


0.0587 


0.8478 


6.6569 


60.6997 


0. 0403 


0.2498 


7 


0.0557 


0.5745 


6.6825 


21.3200 


0.0348 


0.1672 


8 


0.0522 


0.3709 


6.0419 


8.2223 


0.0301 


0.1153 


9 


0.0482 


0.2307 


5.0782 


6.9855 


0.0261 


0.0832 


10 


0.0439 


0 . 1400 


4.0514 


9.3412 


0.0226 


0.0636 


11 


0.0396 


0.0847 


3.1097 


11.6978 


0.0194 


0.0516 


12 


0.0353 


0.0529 


2.3160 


12.9672 


0.0166 


0.0442 


13 


0.0312 


0.0361 


1.6824 


13. 1242 


0.0142 


0.0394 


14 


0.0273 


0.0282 


1.1960 


12.4778 


0.0120 


0. 0360 


15 


0.0237 


0.02 55 


0.8337 


11.3650 


0.0100 


0.0332 


16 


0.0205 


0.0254 


0.5710 


10.0516 


0.0083 


0.0308 


17 


0.0175 


0.0264 


0.3855 


8.7171 


0.0069 


0.0285 


18 


0.0149 


0.0275 


0.2581 


7.4696 


0.0056 


0.0261 


19 


0.0126 


0.0283 


0.1736 


6.3652 


0.0045 


0.0237 


20 


0.0105 


0.0286 


0.1199 


5.4260 


0.0036 


0.0212 


21 


0. 0087 


0.0282 


0.0877 


4.6526 


0.0028 


0.0188 


22 


0.0072 


0.0273 


0.0701 


4.0336 


0.0022 


0.0165 


23 


0.0058 


0.0259 


0.0619 


3.5512 


0.0017 


0.0142 


24 


0.0047 


0.0241 


0.0595 


3.1854 


0.0013 


0.0121 


25 


0.0038 


0.0221 


0.0601 


2.9159 


0.0010 


0.0102 


26 


0.0030 


0.0199 


0.0619 


2.7237 


0.0007 


0.0084 


27 


0.0024 


0.0177 


0.0639 


2.5918 


0.0005 


0.0069 


28 


0.0019 


0.0155 


0.0654 


2.5057 


0.0004 


0.0056 


29 


0.0015 


0.0133 


0.0659 


2.4529 


0.0003 


0.0044 


30 


0.0012 


0.0114 


0.0654 


2.4236 


0 . 0002 


0.0035 


31 


0.0009 


0.0096 


0.0640 


2.4099 


0.0001 


0.0027 


32 


0.0008 


0.0081 


0.0618 


2.4059 


0.0001 


0.0021 


33 


0.0006 


0.0067 


0.0590 


2.4071 


0.0001 


0.0016 


34 


0.0006 


0.0056 


0.0557 


2.4105 


0.0001 


0.0012 


35 


0.0005 


0.0046 


0.0523 


2.4141 


0.0001 


0.0008 


36 


0.0005 


0.0039 


0.0489 


2.4167 


0.0000 


0.0006 


37 


0.0004 


0.0033 


0.0455 


2.4177 


0.0000 


0.0004 


38 


0.0004 


0.0028 


0.0424 


2.4170 


0. 0000 


0.0003 


39 


0.0004 


0.0024 


0.0396 


2.4146 


0.0000 


0.0002 


40 


0.0004 


0.0022 


0.0371 


2.4109 


0.0000 


0.0002 


41 


0.0004 


0.0020 


0.0349 


2.4062 


0.0000 


0.0001 


42 


0.0004 


0.0018 


0.0330 


2.4008 


0.0000 


0.0001 


43 


0.0004 


0.0017 


0.0315 


2.3951 


0.0000 


0.0001 


44 


0.0004 


0.0017 


0.0303 


2.3894 


0.0000 


0.0001 


45 


0.0004 


0.0017 


0.0293 


2.3839 


0.0000 


O.OOOl 


46 


0.0004 


0.0016 


0.0285 


2.3787 


0.0000 


0.0001 


47 


0.0004 


0.0016 


0.0279 


2.3741 


0.0000 


0.0001 


48 


0.0004 


0.0016 


0.0275 


2.3700 


OiOOOO 


0.0001 


49 


0.0004 


0.0016 


0.0272 


2.3664 


0.0000 


0.0000 


50 


0.0004 


0.0016 


0.0270 


2.3635 


0.0000 


0.0000 



PK represents the diagonal elements of the theoretical 
covariance matrix for the optimal filter. 

PR represents the diagonal elements of the theoretical 
covariance matrix for the reduced filter. 



96 



Table V-4 



EXPERIMENTAL VARIANCE OF ESTIMATION ERROR 



K 


VK ( 1 , 1 » 


VK(2,2» 


VK ( 3 » 3 ) 


VK(4,4) 


VR(l,n 


VR(2,2» 


1 


0. 1510 


1,9066 


1.9733 


4.0880 


0.1486 


1.9041 


2 


0,1022 


1. 8632 


1,7310 


225.4341 


0. 1008 


1.9624 


3 


0.0807 


1.7518 


2.2373 


326,4941 


0.0801 


2.0096 


4 


0.0667 


1.5112 


3.8145 


257.3462 


0.0684 


1.9437 


5 


0,0677 


1.1743 


5.6299 


138.0874 


0. 0713 


1.7625 


6 


0.0529 


0.8372 


6.0810 


60,2167 


0,0664 


1.5128 


7 


0.0518 


0.5928 


6.5273 


23.2761 


0.0728 


1.2314 


8 


0.0487 


0.3939 


6.0203 


8.7271 


0. 0783 


0.9501 


9 


0.0464 


0.2513 


5,1600 


7.2510 


0.0825 


0.6993 


10 


0.0457 


0. 1612 


4,4359 


9.6855 


0, 0860 


0,4897 


11 


0.0420 


0.1002 


3.5490 


12.4134 


0,0863 


0, 3264 


12 


0.0391 


0, 0617 


2,7683 


14.0120 


0,0890 


0.2020 


13 


0.0368 


0.0401 


2.0579 


14. 2776 


0, 0877 


0.1218 


14 


0.0297 


0,0300 


1.3831 


14.0637 


0. 0830 


0. 0696 


15 


0.0264 


0. 0256 


0.9971 


11.9529 


0,0804 


0.0388 


16 


0.0226 


0.0247 


0.6811 


12.3539 


0. 0756 


0.0270 


17 


0.0188 


0.0260 


0,4289 


9,7735 


0.0704 


0.0274 


16 


0.0153 


0.0261 


0,2742 


7.4905 


0.0639 


0.0335 


19 


0.0131 


0.0271 


0,1830 


6,0835 


0,0579 


0, 0440 


20 


0.0112 


0,0269 


0.1292 


5.5882 


0.0518 


0,0552 


21 


0.0094 


0.0274 


0,0892 


5.2477 


0.0456 


0.0661 


22 


0.0074 


0.0252 


0,0707 


4, 1038 


0. 0396 


0.0745 


23 


0.0061 


0,0240 


0,0651 


3.5547 


0.0340 


0.0811 


24 


0.0047 


0.0216 


0.0634 


2.8311 


0.0287 


0.0846 


25 


0.0039 


0,0202 


0.0596 


3.2233 


0.0239 


0.0858 


26 


0,0032 


0.0185 


0.0569 


2.5204 


0.0196 


0.0844 


27 


0.0025 


0.0166 


0,0547 


2.7723 


0.0159 


0.0809 


28 


0. 0020 


0.0149 


0.0546 


2.8033 


0.0126 


0.0756 


29 


0,0016 


0.0130 


0,0550 


2.6636 


0. 0098 


0.0688 


30 


0.0013 


0,0114 


0.0581 


2.0651 


0.0075 


0.0612 


31 


0. 0010 


0.0098 


0.0584 


2.4213 


0.0057 


0.0532 


32 


0.0008 


0.0083 


0,0567 


2. 7746 


0, 0042 


0.0453 


33 


0.0007 


0,0071 


0.0539 


2.5252 


0.0031 


0,0379 


34 


0.0006 


0.0060 


0.0490 


2.1855 


0.0022 


0.0313 


35 


0.0005 


0.0050 


0.0435 


2.5705 


0.0016 


0.0254 


36 


0. 0005 


0.0042 


0.0420 


2.2536 


0.0011 


0.0204 


37 


0.0004 


0.0035 


0.0437 


2.5342 


0, 0008 


0.0162 


38 


0,0004 


0.0029 


0.0446 


2.3070 


0,0006 


0.0127 


39 


0.0004 


0.0025 


0.0421 


2.0512 


0,0005 


0.0098 


40 


0.0004 


0,0022 


0,0366 


2.4148 


0,0005 


0,0075 


41 


0.0004 


0.0019 


0.0300 


2.2901 


0,0004 


0.0056 


42 


0,0004 


0.0018 


0,0256 


2.5146 


0.0004 


0.0043 


43 


0,0004 


0.0017 


0.0241 


2,2052 


0.0004 


0.0033 


44 


0.0004 


0.0016 


0.0261 


2.3978 


0,0004 


0.0025 


45 


0.0004 


0.0015 


0.0300 


1.9640 


0.0005 


0.0020 


46 


0.0004 


0.0014 


0.0317 


1,9171 


0.0005 


0.0017 


47 


0,0004 


0,0014 


0.0307 


2.2013 


0.0005 


0.0015 


48 


0,0004 


0,0014 


0.0299 


2.3214 


0. 0005 


0.0015 


49 


0.0004 


0,0015 


0.0281 


2.2047 


0.0005 


0.0015 


50 


0.0004 


0,0016 


0.0229 


2.8952 


0.0005 


0.0016 



VK represents the diagonal elements of the experimental 
covariance matrix for the optimal process. 

VR represents the diagonal elements of the experimental 
covariance matrix for the reduced process. 
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Table V-5 



EXPFRIMENTAL MEANS OF ESTIMATION ERROR 



K 


MEANEl 


MEANE2 


MEANE3 


MEANE4 


MEANEIR 


MEANE2R 


1 


0. 0496 


0.0500 


-0.0350 


-0.1543 


0.0496 


0.0500 


2 


0.0372 


0.0392 


-0.0572 


-0.6276 


0.0349 


-0.1157 


3 


0.0356 


0.0307 


-0.0903 


-0.7992 


0.0266 


-0.2315 


4 


0.0380 


0.0271 


-0. 1340 


-0. 7676 


0.0193 


-0.3045 


5 


0.0395 


0.0203 


-0. 1637 


-0.3713 


0.0100 


-0.3441 


6 


0.0333 


-0.0030 


-0.1219 


-0.0534 


-0. 0047 


-0.3578 


7 


0.0209 


-0.0146 


-0.0909 


0.2345 


-0.0167 


-0.3498 


8 


0.0044 


-0.0565 


0.1352 


0.1478 


-0.0403 


-0.3249 


9 


-0.0017 


-0.0536 


0. 1711 


0.0767 


-0.0515 


-0.2894 


10 


0.0049 


-0.0373 


0.1013 


0. 1589 


-0.0548 


-0.2527 


11 


0.0142 


-0.0265 


0.0168 


-0.0124 


-0.0556 


-0.2173 


12 


0.0021 


-0.0281 


0.0951 


-0. 0444 


-0.0653 


-0.1749 


13 


-0.0028 


-0.0232 


0.1136 


-0.1564 


-0.0704 


-0.1356 


14 


-0.0072 


-0.0168 


0.1219 


-0.3115 


-0.0742 


-0.0982 


15 


-0.0086 


-0.0108 


0. 1078 


-0.4010 


-0.0757 


-0.0651 


16 


-0.0058 


-0.0083 


0.0729 


-0.3752 


-0.0745 


-0.0382 


17 


-0.0028 


-0.0080 


0.0511 


0.0015 


-0.0724 


-0.0154 


18 


0.0020 


-0.0112 


0.0313 


-0.0783 


-0.0691 


0.0028 


19 


0.0055 


-0.0150 


0.0153 


-0.0996 


-0.0657 


0.0183 


20 


0.0044 


-0.0137 


0.0109 


-0.1052 


-0.0635 


0.0337 


21 


0. 0040 


-0.0136 


0.0142 


0.2029 


-0. 0606 


0.0456 


22 


0. 0014 


-0. 0094 


0.0236 


0.1062 


-0.0581 


0.0570 


23 


0.0022 


-0.0103 


0.0273 


0.0610 


-0.0543 


0.0640 


24 


0.0017 


-0.0089 


0.026? 


-0.0824 


-0.0506 


0.0699 


25 


0.0011 


-0.0072 


0.0235 


-0.0285 


-0.0469 


0.0738 


26 


0.0008 


-0.0063 


0.0215 


-0.0504 


-0.0430 


0.0758 


27 


0.0007 


-0.0056 


0.0207 


-0.0033 


-0.0390 


0.0762 


28 


0.0006 


-0.0048 


0.0263 


0. 1783 


-0.0351 


0.0755 


29 


0.0001 


-0.0027 


0.0329 


0.1111 


-0.0314 


0.0744 


30 


0.0002 


-0.0014 


0.0292 


-0.2192 


-0.0277 


0.0723 


31 


0.0003 


-0.0006 


0.0215 


-0.1151 


-0.0241 


0.0693 


32 


0.0004 


0.0002 


0.0170 


-0.0758 


-0.0207 


0.0656 


33 


0.0005 


0.0008 


0.0080 


-0. 2516 


-0.0175 


0.0613 


34 


0.0004 


0.0011 


0.0012 


-0. 0467 


-0.0146 


0.0564 


35 


0.0004 


0.0012 


0.0022 


0.0639 


-0.0119 


0.0513 


36 


0.0005 


0.0014 


0.0016 


-0.0675 


-0.0094 


0.0463 


37 


0.0006 


0.0014 


-0.0011 


-0.0404 


-0.0072 


0.0413 


38 


0.0006 


0.0012 


-0.0069 


-0.1685 


-0.0053 


0.0363 


39 


0.0006 


0.0007 


-0.0093 


0.0404 


-0.0036 


0.0313 


40 


0.0006 


0.0003 


-0.0057 


0.0877 


-0.0022 


0.0266 


41 


0.0007 


0.0002 


-0.0006 


0. 1207 


-0.0009 


0.0225 


42 


0.0007 


0.0003 


0.0027 


0.0151 


0.0001 


0.0139 


43 


0.0008 


0.0004 


0.0038 


0.0374 


0.0010 


0.0157 


44 


0.0008 


0.0006 


0.0026 


-0. 0628 


0.0017 


0.0129 


45 


0.0008 


0.0007 


0.0004 


-0.0407 


0.0022 


0.0103 


46 


0.0008 


0.0007 


0.0035 


0.1341 


0.0027 


0.0080 


47 


0.0007 


0.0010 


0.0075 


0.0162 


0. 0031 


0.0062 


48 


0.0008 


0.0014 


0.0069 


-0.0321 


0. 0033 


0.0048 


49 


0.0009 


0.0017 


0.0066 


0.0174 


0.0035 


0.0036 


50 


0.0010 


0.0020 


0.0053 


-0.0515 


0.0037 


0.0026 



MEANEl through MEANE4 represent the mean of estimation 
error in the complete state vector by means of optimal 
filtering. 

MEANEIR and MEANE2R represent the mean of estimation 
error in states and X 2 by means of a reduced filter, 
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Variance of estimation error 
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Fig. 5.3 - Graphical plot of the variance of estima- 
tion error in x-^ vs. time. Curves a and b show the 
theoretical variance of estimation error for the 
reduced and full filters respectively. Curves c and 
d show the experimental variance for the full and 
reduced processes respectively. 
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Variance of estimation error 
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Fig. 5.4 - Graphical plot of the variance of estima- 
tion error in vs. time. Curves a and b show the 
theoretical variance of estimation error for the 
reduced and full filters respectively. Curves c and 
d show the experimental variance for the full and 
reduced processes respectively. 
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Mean of estimation error 
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Fig. 5.5 - Graphical plot of the means of estimation 
error in X]^ vs. time. Curve a shows the mean of 
estimation error for the optimal Kalman filter, and 
curve b shows that for the reduced filter. 
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ig. 5.6 - Graphical plot of the means of estimation 
error in X 2 vs. time. Curve a shows the mean of 
estimation error for the optimal Kalman filter, and 
curve b shows that for the reduced filter. 
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103 



system and filter for the dominant -root implementation. 



VI. DISCUSSION OF RESULTS 



The three different approaches presented in this 
thesis have as a common objective a decrease in the compu- 
tational requirements involved in handling systems of high 
dimension, based upon the fact that only a preselected re- 
duced number (r) of states is of significant importance in 
the estimation process. 

For the particular example chosen as an illustration 
for the three methods, the partitioning idea discussed in 
Chapter III, although it is the easiest to apply, appears 
to give the least attractive results in the reduced esti- 
mation of states and X 2 where a small bias is produced. 
This seems to be a direct consequence of the degradation 
introduced when (which in this example has relative 

large elements as compared with was not considered to 

be a part of the dynamics of the primary system. 

In the pseudo-inversion approach, implementation of 
four scalar filters, plus the selection of a proper set of 
matrices and its pseudo-inverses was necessary, but 
since it takes less computing capacity to manipulate four 
scalar filters than one (4 x 4) filter, the estimation 
process was clearly simplified. 
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The selection of the partitioning matrices is arbi- 
trary for any particular problem, therefore for this ex- 



and 



with 



4 - 



another selection 


was made as 




^1 


1 1 


-1 -1 


> 








1 






^2 "" 


1 1 








r 








^3 




-1 1 


9 






1 -1 


1 -1 






f 0.25' 




'0.25'i 


fo.251 


1 r 0.251 
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I 0.25! 
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,-0.25 


i-0.25^ ’ 


A 2 = 


0.25 


’ ^3 j-0.25 


’ \ , 0.25 


i-0.25j 




0.25 1 


[0.25 


i-0.2 5 



This particular choice, however, leads to four fil- 
ters for which the sub-optimal variance of estimation er- 
ror, both theoretical and experimental, increases linearly 
with time, which is undesireable . Therefore this parti- 
cular partitioning was not used. 

The plant for all three examples was selected with a 
pair of clearly dominant complex roots in order to obtain 
a good second-order approximation of the complete process 
from which the reduced filter could be implemented using 
the method presented in Chapter V. 
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APPENDIX A 



MATRIX PSEUDO INVERSION 

A review of matrix pseudo-inversion is presented in ■ 
this appendix as an extension to the idea of pseudo-inver- 
sion introduced in Chapter IV. 

A matrix A"*~ is said to be the pseudo-inverse of a rec 
tangular matrix A if 

A A^A = A (A-1) 

where if A is a square non-singular matrix, i.e., A~^ ex- 
ists, then 



A A’^A = A 

is a trivial solution since A A~^ = I. 

Another important definition may be introduced using 
a set of four matrix equations 



A a"*" A 


= A , 


(A-2a) 


< 

< 

< 


= A+ , 


(A-2b) 


A 


= aa'*'. 


(A-2c) 


! 

A+Af 


= A+A 


(A-2d) 
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This definition implies that the inverse of a matrix A is 
established as the unique pseudo-inverse if and only if it 
satisfies the set of Eqs. (A-2). Furthermore, as a conse- 
quence of the latter definition, the pseudo-inverse among 
other properties has the following: 



A = A 






a"*” = a when A is a non-singular matrix. 



and 



+ T 

0 = 0 . 



The pseudo-inverse A"*” of a matrix A may be shown to 
be unique as follows. 

Let A^ be a set of pseudo-inverse of A^ where 
i = 1,2 . . . M, and assume that the A^ are chosen so that 

M 

AtA. = I . (A-3) 

i=l 

Further assume that another set of pseudo-inverse ^ A^ 
exists for each A^, and that this new set of also sat- 
isfies 







Y^A. 
1 1 



i=l 



I . 



(A-4) 
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Then it follows by proper use of Eqs. (A-3) and (A-4) that 



M 



)_' (A+ - YpA, - 



0 



i=l 



and since the choice of A^ is arbitrary, it follows that 

= A^ with i = 1,2 ..., M 

and the uniqueness property is established. 

For more complete treatment of this subject see 
Penrose (Ref. 8), Werther (Ref. 9) and Deutsch (Ref. 10). 
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APPENDIX B 



IMPLEMENTATION OF A REDUCED SUBCPTIMAL LINEAR FILTER 
FCR DISCRETE DYNAMICS PROCESS, BY PARTITIONING THE 
DYNAMICS CF THE ORIGINAL PROCESS. THE DEGRADATION ERR 
CP OF THE FILTER IS COMPUTED AND USED AS A 
MEANS OF CORRECTING THE FILTER VARIANCES. A MONTE CAR 
LC SIMULATICN OF THE PROCESS IS PERFORMED IN CRDER TO 
OBTAIN APROXIMATED VALUES FOR THE MEANS AND COVARIANCE 
CF ESTIMATION ERROR. 
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DIMENSION GAM(4,1) ,GAMR(2t n,PHI(A,A) ,PHIR(2,2I,A(A,4) 

1, H( 1 ,4) ,B( A,l) ,HT(A,1 ) , ATU.A) ,CR<2*2) ,VW( 1,1 ) ,Q(4, A) 

2, PK(4,4) ,HR( l,2),HTR(2, 1) ,G(4,I ),U(4,4» ,PKIR(2,2I 

3, PK1(4,4I, PKR(2,2) ,VWR(1, l),GR( 2, 1 ) , PKOR ( 2 , 2 I ,COR ( 2, 2 ) 

4, E(2,2) . AX 1(4,4) ,AX2( 4,4) ,AX3 (4,4) ,BX1 (4,1), BX2 (4,1) 

5, BX3(4, i),GAMU(4,l) ,DX2( 1,4) ,DEG(2, 2) ,DX3(2,2) ,UR(2,2) 

6, CX4(2,2),A1(2,2) ,A2(2,2) , OX 1 ( 2 , 1 ) , GAMUR ( 2 , 1 ) ,CX3( 2 . 1 ) 

7, CX1( 1, 1),CX2( 1,1) ,R(1,1) ,Y1(50,4) ,Y2(5C,4) ,Y3(50,2) 

8, Y4(50,2),Y5(50,2) ,Y6(50,4),Y7( 50, 2),Y8( 50,4) ,Y9( 50,2) 
S,Y10(50,2) ,Yll( 50,l),Y12(50,l) ,XE(4),XER(2) ,SE1(50) , 
1SE2(5C),SE3(50),SE4(50),SE1R(50),SE2R(50),SV11(50). 
2SV22(50 ),SV33(50),SV44(50),SVRll(50),SVR22{50),X( 50), 
3LL(1),MM(1 ),LLL(1) , MMM ( 1 ) ,P01 ( 50 ) , P02 ( 50 ) , VOl ( 50 ) 
4,VD2(5C),PR1(50),PR2(50),VR1( 50),VR2( 50),ED1(50), 
5EC2(50) ,ER1(50),ER2(50) 

REAL*8 IT( 12)/* VARIANCES OF ESTIMATION ERROR FOR X(l) 

1 LARA BCXll PART. CORR. E=P •/ 

REAL*8 IP( 12)/* VARIANCES OF ESTIMATION ERROR FOR X(2) 

1 LARA BOXll PART. CORR. E=P »/ 

REALMS IE(12)/'MEAN OF ESTIMATION ERROR FOR X(l) 

1 LARA BOXll PART. CORR. E=P •/ 

REAL*8 IF(12)/*MEAN OF ESTIMATION ERROR FCR X(2) 

1 LARA BOXll PART. CORR. E=P •/ 

REAL LAB/4H / 

DATA NT,NS ,NE,M,MR,N,NR,NZ, I X, T/40 , 50 , 2C0 , 1, 1 , 4 , 2 , I , 3 , 
1 C . 0 5 / 

DATA A(1,1),A(1,2),A(1,3),A(1;4),A(2,1),A(2,2),A(2,3), 
2A(2,4),A(3,1),A(3,2),A(3,3),A(3,4),A(4,1),A(4,2),A(4,3 
3),A(4,4)/C.,l.,0.,C.,0.,0.,i.,0.,0.,0.,0., 1., -400. ,-32 
40. ,-118. ,-19./ 

DATA B(1,1),B(2,1),B(3,1),B(4,1 )/0. , 0. , 0. , 400. / 

DATA H(l,l),H(l,2),H(l,3),H(l,4)/l.,0.,fi.,0./ 

DATA HT(1, 1),HT(2,1),HT(3,1),HT(4, 1)/1. ,0. ,0. ,0./ 

DATA VW(1,1),RR,VWR(1,1 ),UK/0.01,0.25,0.01,0.5/ 

DATA SVl,Shl,AWl,KADA,LC0R/a.5,0.1 ,0.0,1,5/ 

DATA Sl,S2,S3,S4/0. 707, 1.41, 1.41,2.0/ 

DATA AM1,AM2,AM3,AM4/0.0,0.0,0.0,0.0/ 

WRITE(6,222) 

DC 31 K=1,NS 

SV33(K)=0.C 

SV44(K)=C.C 

SVR22(K)=C .0 

SE3( K)=0.0 

SE4(K)=C.O 

SE2R(K)=0.C 

SV11(K)=0.0 

SV22(K)=0.0 

SVR11(K)=C.0 

SE1(K)=C.C 

SE2(K)=C.C 

SE1R(K)=C.C 

CONTINUE 

MONTE CARLO SIMLLATIGN 
DC 1 L=1,NE 

CALL GAUSS! IX, SI, AMI, XI) 

CALL GAUSS! IX, S2,AM2,X2) 

CALL GAUSS(IX,S3,AM3,X3) 

CALL GAUSS! IX, S4,AM4,X4) 

X1K=X1 

X2K=X2 
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c 

c 

c 






66 



C 

C 

C 



81 

6 

123 



12A 



142 



)(3K=X3 

X4K=XA 



GENERATING PHI AND GAMRiA 

IFIL.GT.DGG TO 6 
CALL UNITILtN) 

CC A9 1=1, NR 
CC 49 J=1,NR 
LR( I ,J)=U( I, J) 

CALL UNITIPHItN) 

CALL UNITIAXlfM 
CALL SMPYI B,T,GAM,N,M,0) 
CALL SMPYI e,T,EX2,N,M,0) 
CALL SMPYI A,T, AT,N,N,0) 
CC 66 1=1, NT 



XL=I 

RP=1.0/XL 
PP1=1.C /IXL+l.C ) 

CALL GMPRDI AXl , AT, A,N,N,N J 
CALL SMPYI A,RP ,AX1 ,N,N,0) 
CALL GMADDIPHI,AX1,PHI,N,N) 
CALL GMPRDI AT,BX2,B,N,N,M) 
CALL SMPYI B,RP1,BX2,N,M,0) 
CALL GMADDIGAM,6X2,GAM,N,M> 
CCNTINUE 
^.RITEI 6 , 13 C) 

CC 4C 1=1. N 

WRITEI6,12 5MPHIII ,J) , J=1,N) 
WPITEI6, 131) 

WPITEI6,125)GAM 

GENERATING Q 

CALL GMPRD IGAM,VW,BX1 ,N,M,M) 
CALL GMTRA IGAM,DX2,N,M) 

CALL GMPRDIBX1,CX2,C,N,M,N) 
V»PITEI6,132) 

CC 81 1=1, N 

KPITEI6,12 5MQII,J),J=1,N) 
KRITEI6.135) 

^«PITEI6,136) 

CENT INUE 
CC 123 1=1, N 
CC 123 J=1,N 
PK1II,J)=0.0 



CONT INUE 
PK1I3,3)=2. 

PK1I4,4)=A. 

PKII 1,1 )=C.5 
PK1I2,2)=2. 

XEIl )=r. 

XEI2)=0.0 
XEI3)=0.0 
XEI4)=C.O 
XERI l)=O.C 
XERI2)=0.0 
CC 124 1=1, NR 
CO 124 J=1,NR 
PKIRII, J)=G.O 
El I, J)=C.O 
CCNTINUE 
PKIRII, 1)=C,5 
PKIR I2,2) = 2.0 
EI1,1)=C.5 
E I2,2)=2.0 
IFIL.GT.DGG TO 51 
WRITEI6, 134) 

CC 142 1=1, N 

WRITEI6.125) IPKII I ,J) , J=1 ,N) 
WRITEI6.133) 

KFITEI6,125)XE 



110 



ooo ooo o> OOO OOO OOO O o^n 



1 



12 

13 
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WPITE( 6,138) 

CC 143 1=1 fNR 

hPITE(6,125) (E( I, J) f J=lfNR) 

CONTINUE 

CC 3 K=1,NS 

CALL GAUSSnX,SVl,0.0,Vl) 

CALL GAUSS ( IX»SW1*AW1»W1 ) 

IFIK.EQ.DGO TO 4 

X1K = PHI (1, l)*Xl*PHn 1, 2)*X2+PHI ( 1 , 3 ) ♦ X3 + PH I ( 1 , 4 ) ♦ X4 
1+GAM(1 ,1 )» (UK+Wl ) 

X2K=PHI(2, l)*Xl+PHK2t2)*X2 + PHI ( 2 »3 ) ♦ X3* PH I ( 2 1 4 ) *X4 
2+GAM(2f 1 )*(UK+W1) 

X3K=PHI (3,1)*X1+PHI (3,2)»X2+PHn3t3)*X3+PHI(3f4)*X4 
3+GAM<3,l)*(UK+Wl) 

X4K=PHI (4, 1 )*XH-PHI<4, 2)*X2+PHI (4, 3)*X3+PHI(4f4)*X4 
4+GAM(4,l)*(UK+Wl) 

MEASUREMENT 

Z=H( 1, 1)*X1K+V1 

GAIN AND COVARIANCES CF EST, ERROR FOR COMPLETE SYSTEM 

CALL GMPRD(PKl,FTf BXl,N,NtNZ) 

CALL GMPRD(H,BXl,CX2f NZ,N,NZ) 

SC=1.C/(CX2( If D+RR) 

CALL SMPYieXlfSCfGfNfNZfO ) 

CALL GMPRDIGfHfAXlf NfNZfN) 

CALL GMSUBIUf AXlfAXlfNfN) 

CALL GMPRDIAXlfPKlfPKfNfNfN) 

CALL GMTRAIPHI fAX2fNfN) 

CALL GMPRD(PK,AX2»AX1,N,N,N) 

CALL GMPRDIPHI ,AXl,AX2f N.NfN) 

CALL GMADD<AX2fCfPKl,N,N) 

USING KALMAN FILTER EQUATION WITH FULL PLANT 

IFIK. GT.DGO TO 7 

CALL GMPRDIUfXEfBXl ,N,N,M) 

GO TO 8 

CALL GMPRDIPHI fXEfBXlfNf N,M) 

CALL SMPYIGAMfUKfGAMUtNfMfO) 

CALL GMADDI BXlfGAMUf BXl ,N,M) 

CALL GMPRDIGfHfAXlf NfMfN) 

CALL GMSUBIUf AXlf AX2fNfN) 

CALL GMPRD(AX2f BXlf BX2,Nf N,M) 

CALL SMPYIGf ZfBXl f NfMfO ) 

CALL GMADD(0XlfBX2f XEfNfl ) 

E1=X1K-XE(1) 

E2=X2K-XE<2) 

E3 = X3K-XE< 3) 

E-!r=:X4K-XE(4) 

GAIN AND CCVAR. OF ESTIMATION ERROR FOR REDUCED SYSTEM 

IFIL. GT.DGC TC 11 
CC 12 I=lfNR 

CC 12 J=lfNR 
PHIRII , J) = PHH I f J) 

CC 13 I=lfNR 
CC 13 J=lfMR 
GAMRIIf J)=GAM( If J) 

CC 14 I=lfNR 
CC 14 J=lfNZ 
HRIJf I )=H(Jf I) 

FTRI If J) = HT( If J ) 

C RECUCED 

CALL GMPRDIGAMRfVWRfDXl fNRfMRfMP) 
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call GMTRA(GAMR,CX2,NR,FR» 

CALL GMPPD (DX1,CX2,CR,NR,MR,NR» 
nc A A 1=1, NR 
CC AA J=1,NR 
AA FKCR (I , J) = PK( I , J ) 

11 call GMPRD (PKIR ,HTR,nXl,NR,NR,NZ) 

CALL GMPR0(HR,DX1,CX2,NZ,NR,NZ) 

CX 2 ( l,n = l. J/(CX2( 1 ,1)+RR) 

CALL GMPRD(0X1,CX2,GR,NR,NZ,NZ) 

CALL GMPRD(GR,HR,DX2,NR,NZ,NR) 

CALL GMSUB (UR, 0X3, CX3, NR , NR ) 

IF(LC0R.E0.1)G0 TO 1000 

IF(K,GT.UG0 TO A1 

CALL GMPRD ( DX3 ,E ,CCR,NR ,NR,NR) 

I ) in CCMINUE 

^1 CALL GMPPD(DX3,PK1R,PKR,NR,NR,NR) 

CALL GMSUB(FKR,PKOR,DEG,NR,NR) 

CALL GMTRA(PHIR,DXA,NR,NR) 

IF (LCOR.EQ.l )GC TO 2000 

IF(K.E0.1)G0 TO 10 

CALL GMPPD(PHIR,E,A1,NR,NR,NR) 

CALL GMPR0(A1,DXA,A2,NR,NR,NR) 

CALL GMPRD(DX3,A2,CCR,NR,NR,NR» 

II CC 20 1=1, NR 
CC 2C J=1,NR 

PKR( I, J» = PKR( I, J»+CCR( I,J) 

20 E( I , J)=CCR( I , J ) 

OOG CCNTINUE 

CALL GMPRD(PKR,CXA,CX3,NR,NR,NR) 

CALL GMPRD ( PH I R, OX 3, CXA, NR, NR, NR) 

CALL GMADD(DXA,CR,PK1R,NR,NR) 

LSING KALMAN FILTER EQUATION WITH REDUCED PLANT 

IFIK.GT.l )G0 TO 37 

CALL GMPRD(UR,XER,CX3,NR,NR,l ) 

GO TO 17 

37 CALL GMPRD(PHIR,XER,CX1,NR,NR, 1) 

CALL SMPY(GAMR,UK,GAMUR,NR,M,0 ) 

CALL GMADD(0X1,GAMUR,CX3,NR, 1) 

7 CALL GMPRD(GR,HP,DX3,NR, 1,NR> 

CALL GMSUB (UR, CX3, CXA, NR, NR) 

CALL GMPRD(DXA,CX3,CX1,NR,NR, 1) 

CALL SMPY( GR,Z,CX3,NR,NZ,0) 

CALL GMA0D(CX3,DX1 ,XER,NR,1) 

E1R=X1K-XER { 1 ) 

E2R=X2K-XER(2) 

X1 = X1K 
X2=X2K 
X3=X3K 
XA=XAK 

ENSEMBLE MEAN OF PK/K AND PK/K REDUCED, AND GAINS. 

IF(L.GT.1)G0 TO 26 
Y1(K,1 )=PK( 1 ,1 ) 

Y2(K,1)=PK( 2,2) 

Y3(K,1)=PK(3,3) 

Y1(K,3)=PKR(1,1) 

Y8(K,1 ) = G( 1,1) 

Y8(K,2)=G(2, 1) 

Y8(K,3)=G( 3,1) 

Y8(K,A)=G( A, 1) 

Y2(K,3)=PKR(2,2) 

Y9(K, 1)=PK( A,A) 

Y1C(K,1 )=GR(1, 1) 

Y1C(K,2)=GP(2,1) 

REAL VARIANCE OF EST, ERROR FOR FULL AND REDUCED SYST. 

26 SV11(K)=SV11(K)+E1**2 

SV22(K)=SV22(K )+E2**2 
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SV33(K»=SV33(KI+E3**2 
SVA4(K > = SVA4(K I+E4*’*2 
SVR11(K»=SVR11 (KJ+ElR+^2 
^^2 SVP22(K)=SVR22(K» + E2R**2 

r'FAN OF ESTIM. ERROR FOR FULL AND REDUCED PLANT 

SEl (KI = SE1 (Kl + El 
SF2(K»=SE2(K» + E2 
SE3( K)=SE3(K»+E3 
SE4( K> = SEA (K J + E4 

sfir(k»=5;eir(k> + eir 
SE2R(K » = SE2R(K ) + F2R 

SYSTEM DEGRADATION 

IFILCOR.EQ.IIGC TO 300C 
IFIL.GT.UGO TO 27 
Y6(K ,1 > = DEG( 1, 1 ) 

Y6(K,2)=CEG(2,2) 

Y6(K,3»=COR( 1,1) 

Y6(K,4)=C0R(2,2) 

33JJ CONTINUE 
27 IFIL.LT.NE )GO TC 94 
DC 88 K1=1,NS 
Y1(K1,2)=SV11(K1)/NE 
Y2(K1,2)=SV22(K1)/NE 
Y3(K1,2)=SV33(K1 )/NE 
Y9(K1,2)=SV44( KD/NE 
Y4(K1,1 )=SE1(K1 )/NE 
Y5(K1,1 )=SE2(K1 )/NE 
Y4(K1,2)=SF1R(K1)/NE 
Y5(K1,2 )=SE2R(K1 )/NE 
Y1(K1,4)=SVR11( Kl) /NE-Y4(K1 ,2)**2 
Y2(K1,4)=SVR22(K1) /NF-Y5(K1,2)*‘2 
Y11(K1,1)=SE3(K1)/NE 
Y12( K1 ,1) = SE4( KD/NE 
88 X(K1) = K1 

94 CONTINUE 
3 CONTINUE 

1 continue 

DC nil K=1,NS 
PD1(K)=Y1( K,l) 

PP1(K)=Y1(K,3) 

VD1(K)=Y1( K,2) 

VR1(K)=Y1( K,4) 

PC2(K)=Y2(K,1) 

PR2(K)=Y2( K,3) 

VC2<K)=Y2( K, 2) 

VR2(K)=Y2( K,4) 

ED1(K)=Y4( K,l) 

EC2(K)=Y5(K,1) 

ERl (K)=Y4( K,2) 

EP2(K)=Y5( K,2) 
nil CONTINUE 

WRITE(6,222) 

WPITE(6, 16T ) 

ViPITE(6,55 5) 

DC 30C K=1,NS,KADA 

(«PITF( 6 ,lll)K,Ye(K,l),Ye(K, 2 ),Y 8 (K, 3 ),Y 8 (K, 4 ),Y 10 (K,l) 
1,Y1C (K,2) 

3)0 CONTINUE 

V<RITE(6,222 ) 

V.RITE(6 ,161 ) 

WRITE! 6,666) 

CC 77 K=1,NS,KACA 

WRITF(6,111)K,Y1(K,1),Y2(K,1),Y3(K,1),Y9(K,1),YI(K,3) 
ltY2(Kf3i 
77 CONTINUE 

WPITE(6,222) 

WR ITE(6, 162) 

WPITE(6,777) 
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eea 

757 

778 

AO 70 

30 

111 

125 

222 

16C 

16 1 

162 

163 

16A 

130 

555 

666 

777 

99<5 

993 

AAA 

A33 

130 

12 1 

133 

12A 

135 

136 

122 



CC 888 K=1,NS,KA0A 

WRITE(6. 111>K,Y1(K,2) tV2(K,2) »Y3(K,2) ,V9(K,2) .YIIK.A) 
1,Y2(K,A) 

CCNTINUE 

WRITE(6,222) 

V*RITE(6,163I 
I^RITE(6,999I 
CC 757 K=1 ,NS, KAOA 

WBITE(6,111)K,YA(K,1) , Y5 ( K , 1) , Y 1 1 ( K , 1 I , Y1 2 ( K , 1 ) ,Y4( K, 2 
U ,Y5(K,2) 

CCM INUE 

IFdCOR.EQ. DGC TO AOOO 
WRITE(6,222> 

VsRITE(6,16AI 

WRITE(6,998» 

CC 778 K=1,NS,KADA 

WPlTE(6,imK,Y6(K, n ,Y6(K,2) fY6(K, 3) ,Y6(K,A I 
CCNTINUE 
CCNT INUE 

IFILCCR.NE.l IWRITE(6,AAA) 




FORMAT (215 ) 



FORMAT (9X, I2»12F10.A/ ) 
FCRMAT(9X, AFIO.A/) 
FCPMATdHl ) 

pnRMAT#//Qy.inAiK,<: 



F0RMAT(//9X, 
FORMAT (//9X, 
1//I 

FCRMAT(//9X, 
F0RMAT(//9X, 
F0RMAT(//9X, 
FORMAT! ///9X 
FCRMAT(///9X 
1,AX,'G(A,1)» 
FCRMAT(///9X 
2»PK(3,3) 3X 
F0RMAT(///9X 
2* VK( 3,3) »t 3X 
F0RMAT(///9X 
1, AX, 'MEANEA* 



GAINS FCR FULL ANC RECUCED FILTERS*//) 
THEORETICAL VARIANCE OF ESTIMATION ERROR* 




1, AX, 'MhANEA* 

FORM AT (///9X,-rv-,:PA,-u\i.,j.,',*,A,-u»4.,c, 'fTAt-ciifi./' 
1,AX,*E(2,2)»//) 

FCRMAT(//9X, 'CORRECTION DUE TO DEGRADATION APPLIED*) 
cr-ou*T#//«v CORRECTION FOR DEGRADATION* ) 

ILJT / % 



rurvnHi % / / i luii \ j\j 

FCRMAT (//9X, *NO CORRECTION 
FCRMATI9X, *PHI MATRIX*/) 
FORMAT! ///cy- tGAMMA MATRIX*) 



rL.f\rtMI \ r 

F0RMAT(///9X, *GAMMA MATRIX*) 
crouATi ///«v , IMTIALIZATION-XO/-1*/) 

* INITIALIZATION-PO/-1*/) 
*VARIANCE OF FORCING NOISE=n.01 



runriMi \ // / 

FCRMAT(///9X 

FORMAT! ///9X 

FCRMAT!///9X 

FCRMAT! ///9X 

F0RMAT!///9X 

STOP 

END 



ini I 1 — 1 ■ / # 

*VARIANCE OF FORCING NOISE=n.01*) 
*VARIANCE OF MEASUREMENT N0ISE=0.25*) 
*Q MATRIX*/) 
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SUeROUTIKE LMT(U,N) 
CI^'ENSION U(4,4) 

CC 1 1=1, N 
CC 1 J=1,N 
U( I , J)=C,0 

IFU.EQ. I)U( I, J» = 1.9 
1 CCNTINUE 
RETURN 
ENC 



SUBROUTINE GAU SS ( I X , S , AR , V ) 
A=C«C 

DO 50 1=1,12 

CALL RANOU( IX, IY,Y ) 

IX = I Y 
50 A=A+Y 

V= (A-6.C )*S+AM 

RETURN 

END 



SUBROUTINE R ANDU ( I X , lY , YFL ) 

IY=IX^65535 
IFUY) 5,6,6 

5 IY=IY+21A74836A7+1 

6 YFL=IY 

YFL=YFL*.4656613E-S 

RETURN 

ENC 



SUBROUTINE LCC ( I , J , IR , N ,M , MS » 

IX = I 
JX=J 

IF (MS-1) 1C, 20, 30 
10 IRX=N*( JX-D + I X 
GO TO 36 

20 IF(IX-JX) 22,24,24 
22 IRX=IX+( JX*JX-JX)/2 
GC TO 36 

24 IRX=JX+( IX+IX-IX )/2 
GC TO 36 

I R y =6 

IFUX-JX) 36,32,36 
32 IRX=IX 
36 IR=IRX 
RETURN 
END 



SUBROUTINE SMPY ( A , C , R , N ,H , MS ) 
CIMENSICN Ad) ,R( I > 

CALL LOC(N,f',IT,N,M,MS) 

DC 1 1 = 1. IT 
1 R(I)=A( I)*C 
RETURN 
END 



115 



SL8R0UTINE M IN V ( A , N ,D , L t M ) 

CINENSION A( 1) ,L( n,M( 1» 

C = l.f! 

NK=-N 

cc er K=i,N 

^K=NK+N 
L(K)=K 
f^(K)=K 
KK=NK+K 
B IGA=A( KK) 

DC J=K,N 
lZ=N*i J-1 ) 

CC 2C I=K,N 
1J=IZ+I 

10 IF< ABSiBlCAt- ^BS(A(Um 15,20,20 
15 BIGA=A( IJ) 

L (K )=I 
R(K)=J 
20 CCMINUE 
J=L( K) 

IF(J-K) 35,35,25 
25 KI=K-N 

CC 3C 1=1, N 
K I=KI+N 
HCLD=-A( KI ) 

J I=K I-K+J 
fl(KI )=A(JI ) 

30 A( JI ) =HCLD 
35 I=N(K) 

IF(I-K» 45,45,38 
38 JP=N*(I-1) 

CC 4C J=1,N 
JK=NK+ J 
JI=JP+J 
hCLD=-A( JK ) 

A( JK)=A(JI ) 

^0 A( JI ) =HOLD 

45 IF(BIGA) 48,46,48 

46 C=0.0 
RETURN 

48 CC 55 1=1, N 

IF(I-K) 50,55,50 
50 IK=NK+I 

A( IK)=A( IK )/(-BIGA) 

55 CCMINUE 
CC 65 1 = 1, N 
IK=NK+I 
HCLD=A( IK) 

IJ=I-N 
CO 65 J=1,N 
U = I J+N 

IF(I-K) 60,65,60 
60 IF(J-K) 62,65,62 
62 KJ=IJ-I-*-K 

A( IJ ) = HGLD*A<KJ )+A( IJ) 

65 CCMINUE 
KJ=K-N 
DC 75 J = 1,N 
KJ=K J+N 

IF(J-K) 70,75,70 
7C A(KJ)=A(KJ )/BIGA 
75 CCMINUE 
C=D*BIGA 
A( KK )=l,r/BIGA 
80 CCMINUE 
K=N 

IfO K=(K-1) 

IF(K) 150, 150, 105 
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IC5 I=L(K) 

IF(I-K) 12r,12C,i:a 
1C8 JC=N*(K-1) 

JR=N*( I-l) 

DO IK J=1,N 
JK=JO+J 
HOLD=A( JK) 

J I = JR+J 
A( JK )=-A(J I » 

110 A(JI) =HCLC 
120 J=M(K) 

IF(J-K) lUCflOO,125 
125 KI=K-N 

DO 120 1=1, N 
KI=KI+N 
HCLO=A( KI ) 

JI=KI-K+J 
A(KI )=-A(J I ) 

13C A(JI) =HOLD 
GO TO lOG 
150 RETURN 
END 



SUBROUTINE GMTR A ( A . R , N, R J 
DIRENSICN A(l) ,R(1) 

IR=0 

DC 10 1=1, N 

IJ=I-N 

DC 1( J = 1,F' 

IJ=IJ+N 

IR=IR+1 

1C R( IR)=A( IJ ) 

RETURN 

END 



SUBROUTINE GMPRD ( A , B , R , N, M, L) 
DIMENSION A(l) ,6(1) ,R(l) 

IR=0 

IK=-M 

DO 1C K=1,L 

IK=IK+M 

DC 1C J=1,N 

IR=IR+1 

JI=J-N 

IB = IK 

R( IR)=0 

DC 10 1=1, M 

JI=JI+N 

IB=IB+1 

R ( IR)=R( IP) ♦AI J I )*B( IB ) 

RETURN 

END 



SUBROUTINE GMSU E ( A , 8 ,R , N, M ) 
DIMENSION A( 1) ,B( 1) ,R( 1) 
NM=N*M 

DO 10 1=1, NM 
10 R( I )=A( I )-B( I ) 

RETURN 

END 



C 



10 
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SteRCUTINE GMADCIAf B,R,N,M) 
DIMENSION A( 1) ,B< U ,R( 1) 
NR=N*M 

DC 1^ I=ltNM 
10 R( I )=A( I ) + B( I ) 

RETURN 

END 
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APPENDIX C 



I NPLEMeNTAT ION OF A PECUCEO SUP-OPTiMAL FILTER USINC, 
VATRIX PSFUOO INVERSION. FOUR ONE DIMENSIONAL FILT- 
ERS ARE PUN IN PARALLEL IN ORDER TO ESTIMATE ONLY TWO 
STATES OF A FOURTH ORDER SYSTEM, A MONTE CARLO SIMULA 
TICN IS PERFORMED IN ORDER TO OBTAIN APPROXIMATED 
VALUES FOP the MEANS AND EXPERIMENTAL COVARIANCES 
OF ESTIMATICN FPROR for the complete and REDUCED PRO- 
CESSES. 



Cl PENSION PHI ( ,GAM( At 1 ),A(A,A),B(A,n,Q(A,A).Hii,A 

n ,HT( A, 1 > , p ( 1, 1 ) ,G( A, I ) ,PK1 ( A ,A I,PK(A,A),Vw(ltljtU(A,A 

2>tVWH(l,l),AXl(A,A),AX2(A,A>,AX3(A.A),AT(A,A),S(A,c,), 
3S1(A,A),S2 (A,a),hX1(a,1), rX2 (A, 1) » BX3 (A, 1 ) t B XA ( A , 1 ) , B X 
A5(A,ntAIl(A,l),AI2(A,U,AI3(A,l>,AlA(A,n,AYl(A,l),AY 
52 (A, 1) , AY3 (A,l ) ,AYA(a, AY5(A. 1 ),GAMU(a, UtAl(ltAl,A2 
t( 1 ,A) , A3( 1 , A) t AA( 1 , A) ,B11 ( 1 ,A» , B12 (I ,A) ,B13( I ,A > , BIAI! 
7,A),GH( 1,A),DX2(1,A>,CY2( ltl)tCX3( 1,1) ,SE1(5 ) ,SF2(5. ) 
8,SE3(5 ),SFA(5' ),SEIR(S ),SE2R(5 ) , SV 1 1 ( 5. ) , S V22 ( 5 ), 
FSV33(5 ) ,S VAA( 5. ) , SVRl 1 ( 5 ) , S VR22 ( 5 < ) , XAX I S ( 5 ),XE(A), 
ILl (2 ) ,MM(2 ) ,LLL( 1 ) ,MMMI L ) ,Y1 ( P , A ) , Y2 ( 5’'> , A ) , Y 3 ( 5'^ , 3) , 
2YAI5 ,2) ,Y5(bj ,2) ,Y6(5 , 1 ) , Y8 ( 5 • , 6 ) ,Y9 ( 5 ,3)tYU(5M) 

3,PD2(5 ) tVDK 5, ) , VD2( 5 • ) , SDK 5 ) t SD2 ( 5 i ) , VRH 5 ),VP2(5 
A' ),ED1(5= ) ,ED2(5 ),ER1(5 ) , EP 2 ( 5 ) , PDl ( 5* ) , X ( 5 ) 

REAL«8 IT( 12)/' VARI ANC6S OF ESTIMATION EPRCR FOR X(l) 

1 LARA BOX 11 •/ 

RFAL»P IP( 12 )/ • VARI ANCES OF FSTIMATION ERROR FOR X(2) 

1 LAP A BOX 11 •/ 

KEAL>^-8 IF(12)/'MEAN OF ESTIMATICN ERROR FOR X(l) 

1 LARA BOX 11 '/ 

PEAL«8 IF(12)/'M6AN OF ESTIMATION FRROR FOR X(2) 

1 LARA DCX 11 •/ 

REAL LAB/AH / 

DATA NTtNS ,NE,M,MR,NtMR,NZ, IX,T/A.. ,5. , 1 , 1 , A , 1 , 1 , 3 , 

It 5/ 



DATA A(l,l),A(l,2),A(l,3)tA(ltA),A(2,l),A(2,2)tA(2,3)t 
2A(2tA),A(2tl)tA{3t2),A(3t3),A{3,A)tA(A,l),A(A,2),A(A,3 
3) ,A(A,A )/t, . , 1. , ., .tl.t 1. , -A .,-3? 
A .,-118., -19./ 



DATA B(1,1),B(2,1),B(3,1),H( 
DATA H( 1 ,1 ) ,H( 1 ,2) ,H(1 ,3) ,H( 
DATA HT ( 1, 1 ) ,HT( 2, 1 ) ,HT( 3, 1 ) 
DATA VW( 1, 1 ) ,R ( 1 , 1 ) ,VWR( 1 ,1 ) 
DATA UK,SV1,SW1,AW1/ .5, .5, 
DATA V1,V2,V3,VA,AM1,AM2,AM3 

DATA' All K lltAi ( 1,2) ,A1( 1 ,3) 
DATA A2U, 1),A2(1,2),A2( 1,3) 
DATA A3(l,l) ,A?(1,2),A3(1,3) 
DATA AA(1, 1),AA( 1,2),AA(1,3) 
DATA AI1(1,1),AII(2,1),AI1(3 
DATA AI2( 1 ,1) , AI2(2,1 ) ,AI2( 3 
DATA AI 3( 1 , 1 ), AI 3( 2,1 ) ,A! 3( 3 
DATA AI 4(1 ,1 ) , AI4( 2 ,n , AI4( 3 
WRITE( 6,22?) 

DC 31 K=1,NS 
SV33(K)=. .* 

SV44(K)= . 



4,l)/»<», .,'.,4 »/ 

1,4)/1.,- ?./ 

,HT(4,1)/1., ., ., ./ 
,KADA/. >!,'>. 25, w. 1,1/ 

,aAa/*. 7.7, 1,41, 1.41,2 

,A1(1,4)/1., ./ 

,A2(l,4)/v.,l., ., ./ 

, A3 (l,4)/~., .,1., ./ 

,A4(1,4)/.)*^' ,, / ,,!./ 

, 1 ) , A I 1 ( 4 , 1 ) / 1 . , t . , . , 

,1) ,AI2(4,1)/D.,1., ., 
, 1 ) ,AI 3(4, 1)/ , , ,1. , 

,1) ,AI 4(4,1 )/. ., ., ., 



. , 



. / 
. / 
, / 

1 ./ 



SE3(K)=i ,r 
SF4( K) =C .. 

S V 1 1 ( K ) = . 
SV22(K)=: .» 



SVRll(K)=v>. i 
SVR22(K )= . 
SEl ( K)=! 
SE2(K)= 

SEIP (K )=t 
SE2R(K) . 

1 CONTINUE 



MCNTE CAPLC SIMULATION 
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CC 1 L=1,NE 

CALL GAUSS ( IX, VI , AMI , XI ) 

CALL GAUSS! IX, V2,AM2,X2I 
CALL GAUSS ( IX, V3,AM3,X3» 

CALL GAUSS (IX, VA,AMA,XA» 

X1K=X1 

X2K=X2 

X2K=X3 

X4K=XA 

GENERATING PHI AND GAMMA 

lEIL.GT.l »GQ TO 6 
CALL UMIT(U,N» 

UR = 1.( 

CALL UNIT! PHI,N» 

CALL UNIT!AX1,N) 

CALL SMPY! B,T,GAH,N,M*C) 

CALL SMPY! B,T,BX2,N,M,G » 

CALL SMPY! A,T,AT,N,N,T ) 

CC 66 1=1, MT 
XL = I 

RP = l.f /XL 
RP1=1.J/!XL + 1.C ) 

CALL GMPRD! AXl ,AT,A,N,N,N» 

CALL SMPY! A,RP ,AX1,N,N,0» 

CALL GMADD!PHI , AXl ,PHI ,N,NI 
CALL GMPRD!AT,BX2,B,N,N,M» 

CALL SMPY!B,RP1,0X2,N,M, ) 

CALL GMADD!GAM,BX2,GAM,N,M) 

66 CONTINUE 

hRlTE!6,13>. » 

DC AC 1=1, N 

WRITE!6,125MPHI!I,J»,J=1,N» 

WRITE!6,131 » 

WRITEI6 ,125IGAM 

GENERATING Q 

CALL GMPRD!GAM,VW,BX1 ,N,M,M) 

CALL GMTRA!GAM,DX2,N,M ) 

CALL GMPRDIBXl ,DX2,Q,N,M,N) 

WRITE!6,132» 

CC 81 1=1. N 

1 V»RITE!6 ,125M0! I , J» ,J=1,N) 

WRITEI6, 135» 

WR ITE!6,136» 

CONTINUE 
DC 123 1 = 1, N 
DC 123 J=1 ,N 
PKl! I, J»=u. j 
23 continue 

PKl!l,n=f .5 
PK1!2,2)=2. 1 
PK1!3,3 »=2.'< 

PKl!A,A)=A.u 
DC 5A 1=1, N 
CC 54 J=1,N 
54 S( I , JI=PK1 ! I , J » 

XE! 1 )=.,.r 
XC!2)=! .r 
XE! 3)={ 

XE!A»=t .• 

INITIALIZATION VALUES FOR REDUCED FILTERS. 

CALL GMPRDIAl, PKl ,CX2,NR,N,N) 

CALL GMPRD !DX2 ,AI1 ,PKll ,NR,N,NR » 

CALL GMPRD! A2, PKl, DX2, NR, N,N) 

CALL GMPRD !DX2 , AI 2 , PKl 2 , NR , N , NR ) 

CALL GMPRD!A1, XE,YE1,NR,N,NRI 
CALL GMPRD ! A2, XE,YE2,NR,N,NR » 
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CALL GMPRD( A4,PKl,CX2tNRtNf N) 

CALL GMPRD(DX2,AIA,PK1A,NR,N,NR) 

CALL GMPPD( A3,PKl,CX2»NR,NtNl 
CALL GMPRC(DX2,AI3,PK13,NR,N.NR) 

CALL GMPRD( A3, XEf YE3tNR,NfNRl 
CALL GMPPD ( AA, XE,YEA,NR,N,NR) 

TF(L.GT.1)GU to 65 
V»RITE(6,13A) 

DO 1A2 I=l ,N 

142 WRITE(6,125)(PK1( I ,J) , J=1,N) 

WRITE(6,133» 

WRITE(6,125) XE 
65 CCNITINUE 
C 

DC 3 K=1,NS 
C 

call gauss ( IX, SVl,l>.( ,V1 1 
CALL GAUSS(IX,SU1,AW1,W1) 
ifik.eo.dgo to 4 

X1K=PH1 (1, 1 )*X1+PHI (1, 2)*X2+PHI (1, 3)*X3+PHI( 1,4) *X^ 
1+GAM(1,1)*(UK+W1) 

X2K=PHI(2, 1 )*X1+PHI(2,2)*X2+PHI (2,3)*X3+PHI( 2,4)^X4 
2+GAM(2,l)*(UK+Wl) 

X3K = PHI( 3, l)’f'Xl + PHI(3,2)*X2 + PHI ( 3 , 3 ) ♦X3 + PH I ( 3 , 4 ) x-X4 
3+GAM(3, 1)*(UK+W1) 

X4K=PHI (4, 1)#X1 + PHI (4,2)*'X2+PHI (4,3)*X3 + PHI(4,4)^X4 
4+GAMI4, 1)^(UK+W1) 

MEASUREMENT 

4 Z=H( 1, 1 )^X IK+Vl 

GAIN AND COVARIANCES OF EST, ERROR FOR COMPLETE SYSTEM 

CALL GMPRD(FK1,HT,BX1,N,N,NZ) 

CALL GMPRD(H,BX1,CX2,NZ,N,NZ) 

cx2(i, I )=i .*>*/( cx2( 1, 1 )+R( 1, in 

CALL GMPRDIBXl ,CX2 ,G,N,NZ.NZ) 

CALL GMPRD(G,H ,AX1 ,N,NZ,N) 

CALL GMSUB (U, AX1,AX1,N,N) 

CALL GMPRD ( AXl ,PKl ,PK,N,N,N) 

CALL GMTRAIPHI ,AX3,N,N) 

CALL GMPRDIPK, AX3,AX1,N,N,N) 

CALL GMPRDIPHI ,AX1 ,AX2,N,N,N) 

CALL GMADD(AX2,Q,PK1,N,N) 

USING KALMAN FILTER EQUATION WITH FULL PLANT 

IFIK. GT.DGO TO 7 

CALL GMPRD(U,XE,BX1 ,N,N,M) 

GC TO 8 

7 CALL GMPRDIPHI ,XE,BX1,N,N,1) 

CALL SMPYI GAM,UK,GAMU,N,M,0) 

CALL GMADDIBXl ,GAMU,BX1,N,M) 

CALL GMPRD (G,H , AXl ,N, M,N) 

CALL GMSUB (U, AXl, AX2,N,N) 

CALL GMPRD(AX2,eXl,BX2,N,N,M) 

CALL SMPY(G,Z,BX1,N,M,0) 

CALL GMADDIBXl, BX2, XE, N, 1 ) 

E1 = X1K-XEI 1 ) 

E2=X2K-XEI 2) 

E3=X3K-XEI 3 ) 

E4=X4K-XEI 4) 

COMPUTING PHI, GAMMA AND H VALUES FOR REDUCED FILTERS. 

IFIL. GT.l.DR.K.GT. 1 )G0 TO 61 
CALL GMPRDIA1,PHI,611,NR,N,N) 

CALL GMPRDIB11,AI1,PHI1,NR,N,NR) 

CALL GMPRDI A1,GAM,GAM1,NR,N,NR) 

CALL GMPRDIH,AI1,H1,NR,N,NR) 
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CALL GMPRD(A2tPHI»B12,NR,N,N) 

CALL GMPRD( B12tAI2»PHI2,NR,N,NR ) 
CALL GMPRD ( A2,GAM,GAM2 ,NR »NtNR) 
CALL GMPRD(h,AI2»H2»NR,N,NRi 
HT2=H2 

CALL GMPRD (A3,PHI ,B13,NR,N,N) 

CALL GMPRO( 313, AI3, PHI BtNRjN.NR) 
CALL GMPRDI A3,GAM,GAM3,NR,N,NR) 
CALL GMPRD(H,AI 2,H3,MR,N,NR) 
HT3=H3 

CALL GMPRD (A4, Phi ,B14, NR, N,N) 

CALL GMPRD ( B1 4 , A 1 4 , PH I 4 , NR , N , NR ) 
CALL GMPRD ( A4,GAM, GAM4,NR,N,NR ) 
CALL GMPRD(H,AI4,H4,NR,N,NRJ 
PT4=H4 

G REDUCED 

CALL GMPRD (GAMl ,VWR ,CX1 ,NR, MR,MR ) 
CALL GMTRA (GAMl ,CX2 ,NR ,MR > 

CALL GMPRDICXl ,CX2,Q1,NR,MR,NR) 
CALL GMPRD (GAM2 ,VWR ,CXl , NR, MR, MR ) 
CALL GMTRA(GAM2»CX2,NR,MR) 

CALL GMPRD ( CXI ,CX2 ,Q2,NR,MR,NR ) 
CALL GMPRD (GAM3»VWR ,C XI , NR, MR, MR) 
CALL GMTRA(GAM3,CX2,NR,MR ) 

CALL GMPRD(CX1 ,CX2,C3,NR,MR,NR) 
CALL GMPRD (GAM4,VWR ,C X 1 , NR , MR , M R ) 
CALL GMTRA (GAM4,CX2,NR,MR ) 

CALL GMPRDICXI ,CX2 ,C4,NR,MR,NR) 

61 CONTINUE 

SLB-OPTIMAL FILTERS 

CALL GMPRD (PKl 1 ,HT 1,CX1,NR,NR,NZ) 
CALL GMPRD(H1,CX1,CX2,NZ,NR,NZ) 
CX2(1,1) = 1.V(CX2(1,1)+R(1,1)) 
CALL GMPRD(CX1,CX2,G1,NR,NZ,NZ) 
CALL GMPRD(G1,H1,CX1,NR,NZ,NRI 
CALL GMSUB(UR,CX1,CX1,NR,NR) 

CALL GMPRD ( CX 1 , PKl 1 , PR 1 , NR , NR , NR ) 
CALL GMTRAIPHI 1 ,CX2 ,NR ,NR) 

CALL GMPRD(PR1 »CX2,CX1,NR,NR,NR) 
CALL GMPRD (PHI 1 ,CX 1 , CX 2 , NR , NR , NR ) 
CALL GMADD(CX2,Q1,PK11,NR,NR) 
IMK.GT.DGO TO 73 
YP1 = YE1 
YP2=YE2 
YF3=YE3 
YP4=YE4 
GU1=GAM1*UK 
GU2=GAM2«UK 
GU3=GAM3«UK 
GU4=GAM4*UK 
IF(K.EQ.1)GC TO 37 
73 CONTINUE 

CALL GMPRD(AI1,YE1,AY1,N,NR,NR) 
CALL GMPRD(AI2,YE2tAY2»N,NR,NR) 
CALL GMPRD ( AI 3, YE3,AY3,N,NR,NR ) 
CALL GMPRD ( AI4, YE4» AY4,N,NR,NR ) 
CALL GMADD(AY1,AY2,AY1 ,N,NR) 

CALL GMADO( AY3,AY4,AY3,N,NR) 

CALL GMADD(AY3,AY1,AY1,N,NR) 

CALL GMPRD(PHI ,AY1,AY5,N,N,NR) 
CALL GMPRD(A1, AY5,YP1,NR,N,NR) 
CALL GMPRD (A2 , AY5tYP2tNR, N, NR) 
CALL GMPRD(A3,AY5,YP3,NR,N,NR) 
CALL GMPRDIA4, AY5»YP4,NR,N,NR) 

YP 1=YP1+GU1 
YP2=YP2 + GU2 
YF3=YP3+GU3 
YP4=YP4+GU4 
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37 CrMINUe 

CALL GMPRLM A I 1 , YPl , AY 1 , N , NR , NR ) 
CALL GMPRD ( M 2 , YP2 , AY2 ,N,NR ,NR > 
CALL GMPPD( AI3,YP3, AY^,N,NRfNR ) 
CALL oMPRD(AIA,YP-t,AY4,N,NR,NR) 
CALL GMACOC AYl ,AY2tAYl ,N,NR) 

CALL GMAUU ( AY3 , AYA , AY3, N,NR ) 

CALL GMADIH AYl ,AY3,AY2,NfNR) 

FILTLR 1 

CALL GMPKD(G1,H,GH,NR,NR,N) 

CALL GMPRD (GH, AY2, GA,NR,N,NR ) 

GZ=G1*Z 

YE 1 = GZ+YP1 -GA 

FILTER 2 

CALL GMPPD(PK12,HT2,CX1,NR,NR,NZ> 
CALL GMPRD (H2, CXI, CXZjNZ, NR, NZ) 
CX2( 1, 1 ) = 1 . ,/(CX2( 1, 1 ) +R( 1, 1) ) 
CALL GMPRD (CXI , C X2 , G2 , N R , NZ , NZ ) 
CALL GMPRD (G2,H2,C X1,NR,NZ,NR ) 
CALL GMSUB (UR,CX1,CX1,NR,NR> 

CALL GMPRD ( C X 1 , PK I 2 , PR2 , NR , NR , NR ) 
CALL GMTRA(PHI2,CX2,NP,NP ) 

CALL GMPRD(PR2 ,CX2,CX1 , NR, NR, NR ) 
CALL GMPRD (PH I 2, CXI, CX 2, NR, NR, NR) 
CALL (;mADD(CX2,Q2,PK12,NR,NR) 

CALL GMPRD ( G2 , h , GH , NR , NR , N ) 

CALL GMPRD (GH, AY2, GA,NR,N,NR ) 

GZ=G2>?Z 

YE2=GZ+YP2-GA 

FILTER 3 

CALL GMPRD (PK13,HT3,CX1,NK,NR,NZ) 
CALL GMPRD (H3,CX1 ,CX2,NZ,NR,NZ) 
CX2(1,1) = 1. i/(CX2(l,l)+R(l,l)) 
CALL GMPRD(CX1,CX2,G3,NR,NZ,NZ ) 
CALL GMPRD (G3,H3, CXI, NR, NZ,NK) 
CALL GMSUB (HR, CXI , CXI, NR, NR ) 

CALL GMPRD (CXI , PK13 , PR3,NR,NR,NR ) 
CALL GMTRA(PHI3,CX2,NR,NR) 

CALL GMPRD(PR3 ,CX2,CX1,NP,NR,NR ) 
CALL GMPRD ( PH I 3 , CX 1 , CX2 , NR , NR , NR ) 
CALL GMADD(CX2 ,Q3,PK13,NP,NR) 

CALL GMPRD(G3,H,GH,NR,NR,N) 

CALL GMPRD(GH, AY2,GA,NR,N,NR) 

GZ=G3*Z 

YE3=GZ+YP3-GA 

FILTER - 

CALL GMPRD (PKl A,HTA,CX1,NR,NR,NZ) 
CALL GMPRD (HA,CX1,CX2.NZ,NR,NZ ) 
CX2(l,l) = l.,/(CX2(l,l)+R(l,D) 
CALL GMPRD(CX1,CX2,G4,NR,NZ,NZ ) 
CALL GMPRD(G4,HA,CX1,NP,NZ,NR) 
CALL GMSUB (UR, CXI , CXI, NR, NR) 

CALL GMPRD (CXI , PKl A , PR4 , NR , NR , NR ) 
CX2( 1, I )=PHI4 

CALL GMPRD (PR4,CX2,CX1.NR,NR,NR ) 
CALL GMPRD (PHI4,CX 1,CX2,NR,NR,NR ) 
CALL GMADD(CX2 ♦G4,PK14,NR,NR) 

CALL GMPRD(G4,H,GH,NR,NR,N) 

CALL GMPRD (GH, AY2,GA,NR,N,NR) 

GZ=G4»Z 

YE4=GZ+YP4-GA 



123 



SUBOPl IM/6L ESTIMATES 



C 
C 
C 

CALL GMPRD ( AI 1 , YEl ,BX1 ,N,NR,NR ) 

CALL GMPRO( AI2,YE2,BX2,N,NR,NR) 

CALL GMPRH ( AI3 ,YE3t BX4,N,NR,NR ) 

CALL GMPROf AI4,YE4,BX5fN,NR,NR) 

CALL GMAOO( BX4.BX5.BX4,N,NR) 

CALL GMAOD (3X1 ,BX2,BX3,N,NR) 

CALL GMADD(3X4.BX3fBX3fNfNR) 

E1R=X1K-BX3(1, 1 ) 

E2R=X2K-BX3( 2, 1 ) 

X1=XIK 
X2 = X2K 
X3=X3K 
XA=XAK 
C 

C CALCULATING THEORETICAL VARIANCE UF ESTIMATION ERROR 

C 

CALL GMPR0(AI1,G1,BX1,N,NR,NR) 

CALL GMPf<D(AI2 ,G2»BX2,N,NR,NR) 

CALL GMPRD( AI3,G3,BX3,N,NR,NR) 

CALL GMPRD ( AI4,G4, eX4,N,NR,NR ) 

CALL GMADD(BX3,BX4,BX5,N,NR) 

CALL GMArO(BXl,BX2,8X3,N,NR) 

CALL GMAD0(3X3 tBX5, BX3»N,NR) 

CALL GMTRA(BX3,DX2,N,NR) 

CALL GMPRD (BX3»h»AXl,NfNR,N) 

CALL GMSUB (U,AX1 ,AX2,N,M 
CALL GMTRA ( AX2 ,AX1 ,N»N) 

IFIK. EO.l )GO TO 10 

CALL GMPRD(S,AX3fA»N»N,M 
CALL GMPRDIPHI ,A,StN,N,N) 
la CALL GMPRO(AX2,S,A,N»N,N) 

CALL GMPRD(A,AXl,S»NfN,N) 

CALL GMPPO(AX2,Q,A,N»N,N) 

CALL GMPRDI A,AX1 ,S1 ,Nf N,N ) 

CALL GMPRD(BX3,R,BX2,N,NR,NR) 

CALL GMPRD(BX2»DX2,S2,N,NR,N) 

CALL GMADD (S»S1 ,S,N,N) 

CALL GMADD( S,S2fS,N,N) 

C 

C ENSEMBLE MEAN OF PK/K, PK/K SUB-OPTIMAL ANC GAINS 

C 

IFIL. GT.l )GO TO 26 
Y1(K,1)=PK(1 ,1 ) 

Y2(K,1)=PK(2»2) 

Y3 (K,l )=PK (3,3 ) 

Y<?(K,1)=PK(4,4) 

Ye(K,l)=G( 1,1) 

Y8(K,2)=G(2,1) 

Y8(K,3)=G( 3,1) 

YB(K,4)=G(4,1) 

Y8(K,5)=G1 
Y8(K,6) =G2 
Y1(K,4)=S( 1,1 ) 

Y2(K,4)=S(2,2) 

Y2(K,2 )=S( 3,3) 

Y<?(K,2) = S(4,4) 

C 

C REAL COV. OF EST. ERROR FOR FULL AND REDUCED SYSTEMS 

C 

26 SV11(K)=SV11 (K)+E1**2 

SV22(K )=SV22(K )+E2**2 
SV33(K)=SV33(K)+E3**2 
SVAA(K) =SV44(K )+E4*«2 
SVRIKK )=SVR11 (K)+E1P«*2 
SVR22(K)=SVR22 (K)+E2R**2 
C 

C MEAN OF ESTIMATION ERROR FOR FULL AND REDUCED PLANTS. 

C 

SE1(K)=SE1 (K)+E1 
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dS 

55 

3 

1 



3 



77 



383 



757 



C 

C 

C 



35 



I .3 



SF2(K)=St2 (K)+F2 
S63( K) =SE3( K)+E3 
Sf A(K)=SEA(K)+E3 
SE1R(K»=SE13(K»+E1« 

SE2P(K)=SE?R(K)*E2R 

IFd.LT.NE )G0 TO 55 

DC 88 Kl=l ,NS 

Y1(K1,3)=SV] 1( KD/NE 

Y2 (Kl, 3 » = SV22( K1 )/ME 

Y3(K1,3)=SV33( Kll/NE 

Y9 (K If 3 )=S VhA( K l ) /NE 

YA ( Kl , 1 )=S El (Kl l/NE 

Y5(K1,1)=SE2(K1J/NE 

Y6(K1,1)=SE3(K1)/NE 

Y1 1 ( Kl , 1 ) = SEA( Kl )/NE 

YA(K1,2)=SE1R( Kl) /NE 

Y5(K1,2)=SE2R(K1)/NE 

Y1(K1,2)=SVR11(K1) /NF-YA( Kl ,2)’’‘*2 

Y? (K 1, 2 )=SVK22( Kl ) /NE- Y5 ( K 1 , 2 ) “-• 2 

XAXIS( Kl )=K1 

CCNTINUE 

CCNTINUE 

CCNTINUE 

NRITE( 6 f222) 

UR ITE( 6, 16' ) 

UPITE(6,555) 

DC 3'*, K = l,NSfKADA 

UPITE(6,lll)K,Y9(K,l),Y8(K,2»fY8(K,3),Ye(K,A),ytt(K,5), 

1Y8(K,6) 

URITE(6,222) 

URITE(6, 161 ) 

URITE( 6,666) 

CC 77 K=1,NS,KACA 

URITE(6,111)K,Y1(K,1),Y2(K,1 ) ,Y3(K, 1 ) , Y9( K, 1 ) 

UR1TE( 6,222) 

UR1TE(6,162) 

UP1TE(6,777) 

DO 888 K = 1 ,NS, KADA 

UPITE(6,111)K,Y1(K,3),Y2(K,3),Y3(K,3),Y9(K,3),Y1(K,2), 

1Y2(K,2) 

UPITE(6,222) 

WRITE(6, 163) 

UPITE(6,999) 

CO 757 K=1,NS,KADA 

UPITE(6,lll)K,YA(K,l),Y5(K,i),Y6(K,l),Yll(K,l),YA(K,2) 
1 ,Y5( K,2i 
UPITE(6,222) 

URITE(6, 16A) 

WRITE! 6, 525) 

CC 535 K=1 ,NS, KADA 

URITE(6,111)K,Y1(K,A),Y2(K,A),Y3(K,2),Y9(K,2) 

GRAPHICAL OUTPUT 

DC !•' ( K=1,NS 
PD1(K)=Y1(K,1) 

PD2 ( K) =Y2( K, 1 ) 

VCKK ) = Y1( K,3) 

VC2(K)=Y2( K,3) 

SCK K) =Y1( K,A) 

SC2( K ) = Y2( K , A) 

VPl (K)=Y1( K,2) 

VP2( K) =Y2( K,2) 

EC1(K)=YA( K, 1) 

ED2(K)=Y5( K,l) 

EP1(K)=YA( K,2) 

ER2 (K)=Y5( K, 2) 

X(K)=XAXIS(K) 

CONTINUE 

CALL DR AW( 5 ; ,X ,PD1 , I,'" ,LAB, IT , 

CALL DRAU( 5 * ,X ,VD1 ,2,. , LAB, IT, 

CALL DPAW(5 , X , SD 1 , 2 , ? , LAS , I T , 



, , , , ,2,6,12, ,L) 
, , ,.‘,‘,2,6,12, ,L) 
, , , , ,2,6,12, ,L) 
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Ill 

125 

222 

525 

555 

66 6 

777 

999 

13 ' 
131 

132 
13^ 

133 
135 
13o 

16- 
16 1 

162 

163 

16A 



1 



C 



5U 



C 



5 

6 



CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 




»X 

• X 

f X 

• X 

»x . _ . _ 

'.XAXIS, , 

,XAXIS,ER1 ,3,‘, 
,XAXI S.ED2, It 
,XAXIStER2t 3, 



tLABt lEt 't', tOt0fJt2f6t8t‘.'tL » 
tLABtIF t ' f.'ti,ft''j0t2t6t8t'^tLi 
tLABt IFtL't0t0tC't--f2t6t8t 



)tU 



DRAW( 5'^ 

DRAW( 5. 

DRAW( 5; 

DRAW( 5 
DRAW( 5 j 
DRAW( 5J 
DRAW( 53 
DRAU( 5 } 

DRAW(5.* 

F0RMAT(2I5 ) 

FCPMAT(9Xt I2tl2Fl >.A/) 

FCPRAT (9Xt AFin .4/ ) 

FCRMAK IHl ) 

FORMAT (///9Xt • K' t4Xt ' SR( 1 1 1 ) • t 3Xt • SR( 2 t2) • t 3Xt 
l'SR(3t3)'t3Xt' SR(4t4)*//) 

FORMAT (//9X,'K*t5Xt’G( I t 1) ' t 4X t • G( 2 t 1 ) * t 4X t * G ( 3 t 
It *G( ^t 1 » • t 3Xt *G1 ( 1 1 1) • t 3X. *G2( 1 1 1) •//) 

FORMATC //QX t' K* t4X t * PK ( 1 1 1 ) ' t 3X t ' PK (2 1 2 ) * t 3X t ' PK ( 3t 3 P 
lt3Xt 'PK(4t 4)'//) 

FCRMAT(//9Xt‘ ■ ‘ 

ItBXt *VK(4,A) 

FORMAT (///9X 

1 t4Xt ‘MEANE4* 

Ff,RMAT( 9Xt 'PHI MATRIX'/) 

FORMAT (///9Xt ' gamma MATRIX'/) 

FORMAT( ///9Xt' Q MATRIX'/) 

FORMAT ( ///9Xt ' I MT I AL I Z AT I ON-PO /-I ' / ) 

F{ PMAT(///9Xt* IMT I ALIZATION-XO /-!•/) 

rrr. OMAT////OV axiAOTAMrc: f'C cno^Tar- Aini 




1 ) ' t4X 



3Xt *VK(3t3) ' 
'MEANE3' 



FORMAT! ///9Xj ' VARI ANCE “CF ' FORCl NG 'NO I SE=0 .0 1 ' . 

“ 'VARIANCE OF MEASUREMENT N0ISE=1.25') 



) 



FORMAT ( ///9X, . 

FORMAT (//9Xt *GA INS FOR FULL AND REDUCED FILTERS'//) 
FORMAT! //9Xt ' THEORETICAL VARIANCE OF ESTIMATION ERROR* 
1 t IXt FULL PLANT' //) 

FORMAT! //9Xt 'REAL VARIANCE OF ESTIMATION ERROR'//) 
F0RMAT!//9Xt 'MEANS OF ESTIMATION ERROR'//) 

F0PMAT!//9Xt 'SUBOPTIMAL VARIANCE OF ESTIMATION ERROR* 
// ) 



1 //) 
STOP 
END 



SUBROUTINE UNIT!UtN) 
DIMENSICN U!4t4) 

DC 1 I=ltN 
DC 1 J=ltN 
U! It J)=C .0 

IF!J.FQ.I)U!It J)=1.C 

CENT IMUE 

RETURN 

END 



SUBROUTINE GAU SS ! I X t S t AM , V ) 

A=i. *‘ 

CO S'- I = ltl2 

CALL RANOU! IXt lYtY) 

IX=IY 
A = A+Y 

V=!A-6.C)*S+AM 

RETURN 

END 



SLBKOUTINE RAN DU ! I X t I Y t YFL ) 

IY=I X*65539 
IF! lY) 5t6t 6 
I Y=IY*2147483647+1 
YFL=IY 

YFL=YFL*.4656613E-9 

RETURN 
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^3 CC 55 1=1, N 

IF( I-K» 5: ,5 5, 5 . 

5 - IK = NK+ I 

/^( IK) = A( IK)/(-HlGA) 

55 CCMINUE 
CC 65 1=1, N 
I K=NK+ I 
HC LD=A( IK) 

IJ=I-N 

CC 65 J=1,N 

IJ=IJ+N 

IF( I-K ) f ,65,0" 

6 IF(J-K) 62,65,62 
62 KJ=IJ-I+K 

A( IJ ) = HGLD*A(KJ ) + A( IJ ) 
65 CCMINUF 
K J=K-N 
CC 75 J=1,N 
K J=K J+N 

IF(J-K) 7. ,75,7. 

7 A(KJ )=A(KJ )/SIGA 
75 CCMIMUE 

D=C*BIGA 
A(KK )=!.> /BIGA 
CCMIMJE 
K = N 

1. K=(K-1 ) 

IF<K) 15- ,15: , n5 
15 I=L(K) 

IF(I-K ) 12* , 12 , 1 "8 

I 8 JC=N*(K-1) 

JP=N*( I-l) 

CO IK J=1,N 

JK=JC+J 

HCLD=A(JK) 

JI=JR+J 
AUK ) = -A(J I ) 

II A(JI) =HCLD 
12 J=M(K) 

IF(J-K) K ,10» , 12 5 
125 KI=K-N 

CO 130 I=1,N 
KI=KI+N 
HCLD=A( KI) 

JI=KI-K+J 
A(KI ) = -AU I ) 

13 j A( JI) =HCLC 
GO TO li.t 
15' RETURN 
END 



SUBROUTINE GMTRA ( A , R , N , M ) 
DIRENSICN Am ,R(1 ) 

IP=f 

CC ir 1 = 1 , N 

IJ=I-N 

CC 1' j=i,y 

IJ=IJ+N 

IR=IP+1 

1 R(IP) = A(IJ) 

RETURN 

END 



SUBROUTINE GMPRD ( A . 6 , R , N , M, U ) 
DINFNSICN A(l) ,B(l ) ,R(1 ) 

IR=( 

IK=-P 

CO 1< K=1,L 
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SUBRCUTINE LCC ( 1 1 J f IR t N t M , MS I 

IX=I 
JX = J 

IF(MS-l) lv,2'J,3C 
1‘ IRX=N*( JX-1 )+IX 
GC TO 36 

2' IF(IX-JX) 22 , 2 ^, 2 ^ 

22 IPX = IX+( JX*JX-JX )/2 
GC TO 36 

2A IRX=JX+( IX*IX-IX)/2 
GC TG 36 
IPX=0 

IF(IX-JX) 36,32,36 
32 IPX=IX 
36 IP=IRX 
RETURN 
EMC 



SteROUTINE SMPY( A,C,R, N,M,MS ) 
CIMENSIGN A(l) ,R(1 ) 

CALL LOC(N,M,IT,N,M,MS) 

CC 1 1=1, IT 
1 R(I)=A(IKC 
RETURN 
END 



SUBROUTINE M INV ( A, N , D, L , M ) 

DIMENSION A(l) ,L(1) ,M( n 

D=1.0 

NK=-N 

CC 8C K=1,N 

NK=NK+N 

L(K)=K 

M(K) =K 

KK=NK+K 

eiGA=A(KK) 

CG 2t J=K,N 
IZ=N=« ( J-1) 

CO 2"* I=K,N 
IJ=IZ+I 

10 IF( ABS(BIGA)- ABSIAdJM) 15,20,20 
15 BIGA=A(IJ) 

L(K)=I 
M(K)=J 
2-^' CONTINUE 
J=L( K) 

IF(J-K) 35,35,25 
25 KI=K-N 

DC 3f I=1,N 
KI=KI+N 
FCLD=-A(KI ) 

JI=KI“K+J 
A(KI )=A( JI ) 

3} A(JI) =HOLD 
35 I=M(K) 

IF(I-K) 45,45,38 
38 JP=N*(I-1) 

CC 4r J = 1,N 
JK=NK+J 
JI=JP+J 
HCLD=-A(JK ) 

A(JK) = A( JI ) 

4' A(jn =HOLD 

45 IF(BIGA) 48,46,48 

46 D=C.C 
RETURN 
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IK= IK + M 
CC 1» J = 1,N 
IP=IR+1 
J I = J-N 
I0 = IK 
R( IR » = 

CC 1' 1 = 1, R 

JI=JI-*-N 

IR=IB+1 

1 R ( IR »=R ( IR »+A( J I )>^H( IB » 
RETURN 
END 



SUBROUTINE GMSUB ( A , B ,R, N, M ) 
DIRENSICN A(l) ,B(1) ,R(1» 
NR=N*M 

CC If 1=1, NM 
1 R( I )=A( I )-B( I) 

RETURN 

END 



SUBROUTINE GMADO ( A , B , R , N , M ) 
Cl PENSION A(l),R(l),Rm 
NN=N*R 

CO 1C 1 = 1, NM 
1 R(n=Am+E(i) 

RETURN 

ENC 
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